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SUMMARY 


The  following  report,  b<ised  on  the  Ph.D.  dissertation  of  the  first  author 
[141],  presents  the  results  of  an  improved  quasi-static  numerical  simulation  algo¬ 
rithm  developed  to  study  both  mechanical  and  scalar  transport  properties  of  three- 
dimensional  idealized  granular  assemblages  simultaneously.  In  addition,  the  results 
of  an  experimental  investigation  of  these  properties  axe  also  presented  and  compared 
against  the  numerical  predictions.  The  simulation  algorithm  includes  several  new 
techniques,  including  a  shuffling  algorithm  for  the  generation  of  an  initial  random 
packing  of  a  granular  assemblage  and  an  improved  microcell- adjacency  method  to 
accelerate  particle-contact  search.  Furthermore,  a  relaxation  method  is  employed  to 
overcome  a  singularity  in  the  quasi-linear  system  of  equilibrium  equations. 

With  the  objective  of  correlating  scalar  transport  properties  such  as  electri¬ 
cal  conductivity  with  the  mechanical  behavior  of  granular  media,  we  treat  the  granular 
assemblage  as  a  resistor  network,  with  particle  centers  being  nodes  and  interparticle 
contacts  being  resistors,  for  the  purpose  of  computing  the  conductivity. 

The  Reynolds  dilatancy  for  randomly  dense-packed  granular  assemblages  is 
found  to  depend  on  the  interparticle  friction,  at  odds  with  Reynolds’  original  hy¬ 
pothesis.  The  use  of  linear  contact  mechanics  is  found  to  be  valid  neair  the  ideal 
rigid-particle  limit.  Also,  a  strong  correlation  is  found  between  electrical  conductiv¬ 
ity,  stress  and  fabric  tensors,  indicating  that  the  scalar  transport  properties  can  serve 
as  a  useful  macroscopic  probe  for  the  particle-contact  topology  in  granular  media. 

Triajcial  compression  tests,  employing  steel  balls  m  electrically  conductive 
granular  particles  serve  to  confirm  our  simulation  of  both  the  mechanical  and  scalar 
transport  properties,  provided  that  the  electrical  conductivity  calculations  are  based 
on  the  experimental  load-resistance  charaw:teristics  of  individual  contacts.  The  mea¬ 
sured  contact  resistance  between  steel  balls  is  found  to  be  much  higher  than  theoretical 
predictions  based  on  Hertzian  contact,  and  exhibits  a  much  stronger  dependence  on 
normal  load,  possibly  due  to  asperities  and  oxide  films  on  the  steel-ball  surfaces. 
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Chapter  1 


Introduction 


Granular  media  are  materials  composed  of  distinct  particles  which  can  move 
independently  of  one  another  and  which  interact  only  at  highly  localized  interparticle 
contact  regions.  In  general,  a  test  on  real  granular  media  such  as  sand  is  difficult  to 
interpret  since  the  stress  inside  the  sample  can  not  be  measured  directly  and  must  be 
estimated  from  the  boundary  conditions,  although  measurements  of  strain  have  been 
made  possible  by  techniques  such  as  X-ray  photography  technique.  Also,  Dantu  and 
Wakabayjishi  (1957)  suggested  the  use  of  an  photoelastic  material  for  rods  or  discs  in 
order  to  determine  stresses  in  granular  media.  Analysis  of  the  force  distribution  in 
such  a  test  was  first  described  by  De  Josselin  De  Jong  and  Verruijt  (1969),  and  the 
technique  has  been  adopted  by  many  researchers  [45,78,79,102,103,124].  Although 
testing  of  assemblages  of  photoelastic  discs  allows  for  an  accurate  determination  of 
contact  forces,  displacements  and  rotations  of  the  individual  discs,  the  analysis  is 
time  consuming.  Moreover,  the  technique  is  not  as  yet  applicable  to  3-dimensional 
samples. 

While  physical  models  are  certainly  the  ultimate  test  of  any  physical  the¬ 
ory,  numerical  simulation  has  the  advantages  over  real  experiments  in  that  any  mi¬ 
croscopic  information  essential  to  the  understanding  of  the  macroscopic  behavior  of 
these  systems  is  accessible  at  any  stage  of  a  test,  and  “experiments”  can  be  performed 
numerically  that  would  be  very  difficult  physically.  Many  reported  works  show  that 
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numerical  techniques  are  capable  of  reproducing  qualitatively  the  overall  continuum 
mechanical  behavior  of  real  granular  materials  [39,40,126,128,29].  Compared  with 
real  granular  media  such  as  soils,  however,  current  numerical  techniques  are  able  to 
simulate  only  idealized  particle  shapes  such  as  disks,  spheres,  ellipsoids  etc.,  in  a  hm- 
ited  sample  size,  though  the  size  effect  is  partly  overcome  by  the  use  of  periodic  cell 
models. 

Currently,  there  are  mainly  two  classes  of  numerical  technique  employed 
to  simulate  the  quasi-static  mechanics  of  granular  materials,  namely,  dynamic  and 
quasi-static.  The  dynamical  simulation  technique,  often  referred  to  as  the  Distinct 
Element  Method(DEM)  in  the  older  literature,  was  first  developed  by  Cundall  and 
Strack  [39]  and  has  been  widely  employed  since  [40,126,41,128,18,33,10].  However, 
various  artificial  damping  procedures  have  to  be  used  to  suppress  parasitic  particle 
vibrations  in  order  to  achieve  quasi-static  conditions.  Moreover,  it  has  been  noticed 
recently  that  the  algorithm  is  only  conditionally  stable  [10,30].  For  this  reason  among 
others,  a  direct  quasi-static  simulation  has  been  receiving  increased  attention  in  recent 
years  [115,76,26,14,58]. 

Reynolds  dilatancy,  one  of  the  most  fundamental  characteristics  of  granu¬ 
lar  materials,  has  been  accounted  for  in  the  mechanical  modeling  of  granular  flow 
[Reynolds(1895),  Rowe(1962),  Oda(1974a),  Nemat-Nasser(1980)  and  Goddard  et  al. 
(1990b)].  The  factors  influencing  dilatancy  have  been  studied  by  many  investigators 
[120,94,39,103,18,33,14].  The  effects  of  interparticle  sliding  friction  fi  and  of  factors 
such  as  initial  void  ratio  and  state  of  packing  have  been  explored  sporadically,  often 
with  conflicting  conclusions.  As  for  randomly  dense-packed  granular  assemblages, 
Reynolds  suggested  that  friction  should  have  no  effect  other  than  to  stablize  oth¬ 
erwise  unstable  granular  configuration  [107].  Skinner  showed  experimentally  that 
friction  has  little  effect  on  dilatancy  of  random  assemblages  of  spherical  particles 
[120].  On  the  other  hand,  this  effect  has  also  been  investigated  by  means  of  certain 
computer  simulation,  mostly  on  the  random  dense-packed  granular  arrays,  which  in¬ 
terestingly  led  to  the  opposite  conclusion.  In  simulating  the  physical  experiments 
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reported  by  Oda  and  Konishi  [98],  Cundall  et  al.  found  that  dilatancy  depends  on 
interparticle  friction,  the  assembly  with  higher  friction  coefficients  dilated  more  and 
at  a  greater  rate  [39].  Similar  results  have  been  reported  by  Bathurst  et  al.  [18],  and 
Chen  [33]. 

Due  to  their  discrete  nature,  the  behavior  of  the  granular  media  generally 
depends  on  a  variety  of  factors,  such  as  void  ratio,  interparticle  friction,  particle  shape, 
and  microstructural  arrangement  or  “fabric”,  to  name  only  a  few.  Granular  fabric  is 
believed  to  be  one  of  the  most  important  factors  determining  the  overall  mechanical 
response  of  a  medium  to  the  deformation.  Oda  and  Konishi  [94,95,96,98]  performed 
direct  measurement  of  fabric  in  sand  specimens  and  made  many  important  discoveries 
on  the  deformation  mechanisms  of  granular  materials.  However,  such  measurements 
are  difficult  or  tedious  to  perform  experimentally.  It  would,  therefore,  be  highly 
desirable  if  the  granular  fabric  could  be  related  to  and  measured  indirectly  by  means 
of  macroscopic  quantities.  Dynamic  shear  modulus  and  even  the  complete  set  of 
elastic  moduli,  inferred  from  wave  speed  measurements,  has  been  found  to  contain 
direct  information  about  the  internal  fabric  [31,32,66,1,2]. 

On  the  other  hand,  scalar  transport  processes  such  as  electrical  or  thermal 
conduction  through  granular  materials  can  provide  another  such  macroscopic  quan¬ 
tity,  since  the  effective  conductivity  of  granular  materials  depends  not  only  on  the 
conductivities  of  solid  grains  and  interstitial  or  pore  fluid,  but  also  on  the  volume 
fraction  of  solid  particles  (void  ratio)  and  particle  arrangement  or  fabric.  In  fact, 
the  evolution  of  mechanical  anisotropy  of  water  saturated  sands  and  clays  has  been 
studied  in  triaxial  compression  tests  by  monitoring  the  radial  and  axial  electrical 
conductivity  [85,8,4,5].  However,  since  sand  grains  are  themselves  not  electrically 
conductive,  the  current  is  conducted  only  through  the  pore  water.  Therefore,  the 
measured  anisotropy  of  conductivity  mainly  reflects  the  anisotropy  of  the  void  space. 
While  void  space  is  part  of  the  internal  structure,  it  does  not  serve  as  a  good  indicator 
of  granular  contact  topology.  In  particular,  granular  chain  structures,  which  bear  the 
major  load,  and  the  variation  of  interparticle  contact  forces  are  not  fully  captured  by 
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the  conductivity  measurements  of  many  earlier  studies 

Based  on  above  considerations,  the  present  work  is  concerned  with  a  sys¬ 
tem  consisting  of  electrically  conductive  particles  and  an  electrically  nonconductive 
interstitial  fluid.  One  objective  is  to  find  correlation  between  the  scalar  property  and 
mechanical  properties  during  deformation.  The  investigation  includes  both  numeri¬ 
cal  simulation  and  physical  experiments.  The  numerical  simulations  aUow  access  to 
detailed  microstructural  information,  such  as  internal  fabric,  coordination  number, 
local  contact  force  etc..  The  triaxial  compression  tests,  which  employ  steel  balls  as 
conductive  granular  particles  serve  to  validate  the  computer  simulations. 

In  Chapter  2,  some  important  aspects  of  the  theoretical  development  and 
its  application  to  current  investigation  will  be  reviewed  briefly,  including:  (1)  the 
fabric  tensor,  (2)  particle  contact  mechanics  and  nonlinear  elasticity  of  granular  me¬ 
dia,  and  (3)  scalar  transport  through  granular  media.  Chapter  3  provides  a  detailed 
description  of  the  quasi-static  simulation  and  the  various  newly  developed  simulation 
techniques,  while  Chapter  4  covers  the  experimental  aspects  of  the  current  inves¬ 
tigation.  Numerical  simulations,  mainly  aimed  at  the  study  of  the  microstructural 
properties  of  the  media,  are  explored  in  Chapter  5.  Next,  the  results  of  computer 
simulation  on  scalar  transport  through  idealized  granular  assemblages  are  compared 
with  experimental  observations  in  Chapter  6.  Finally,  Chapter  7  summarizes  the 
major  conclusions  of  the  present  study  and  suggestions  for  future  work. 


Chapter  2 


Literature  Review 


2.1  Microstructnre  of  Granular  Media 

It  is  well  accepted  nowadays  that  porosity  or  solid  volume  fraction  alone  is 
not  sufficient  to  characterize  the  geometry  of  the  local  microstructure  of  a  granular 
material,  given  that  two  specimens  of  a  granular  material  such  as  soil,  with  identi¬ 
cal  porosity,  may  possess  quite  different  microstructure  and  behave  mechanically  in 
entirely  different  ways.  In  order  to  understand  the  dependence  of  the  stress-strain  re¬ 
lation  on  microstructure,  additional  geometric  measures  of  local  structure  such  as  the 
geometric  fabric  tensor,  have  been  proposed  by  many  investigators  in  different  fields, 
including  granular  materials,  soil  and  rock  mechanics  [99,100,113,73],  cancellous  or 
spongy  bone  mechanics  [62],  composite  micromechanics  [52]  etc.. 

Oda  (1978)  [99]  introduced  the  concept  of  a  fabric  ellipsoid,  an  ellipsoid 
determined  by  the  three  dimensional  distribution  of  the  unit  normal  to  the  tangential 
contact  planes.  Oda,  Konishi  and  Nemat-Nasser  (1980)  [100]  developed  the  idea  of 
the  fabric  ellipsoid,  equivalent  to  a  second  rank  synunetric  tensor,  and  argued  that, 
after  porosity,  it  is  the  second  best  measure  of  microstructure  in  granular  materials, 
which  appears  to  be  a  matter  of  general  agreement  now.  Following  the  work  of  Oda 
et  al.,  these  second  rank  tensors  are  generally  called  fabric  tensors. 

According  to  Oda  (1978),  fabric  represents  the  spatial  arrangement  of  parti- 
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cles  and  associated  voids.  This  may  includes:  (a)  orientation  fabric,  which  relates  to 
the  inclination  of  a  characteristic  dimension  of  individual  particles  relative  to  a  refer¬ 
ence  direction;  (b)  packing  or  mutual  relation  of  particles,  defined  by  the  probability 
density  function,  E{n),  of  contact  normals  n,  and  the  average  coordination  number 
(the  number  of  contacts  per  particle). 

The  anisotropy  of  granular  materials,  measured  by  the  fabric  tensor,  has 
been  divided  into  two  categories:  (a)  inherent  anisotropy,  a  physical  characteristic 
inherent  in  the  virgin  materials  and  entirely  independent  of  the  applied  strain;  (b) 
induced  anisotropy,  due  exclusively  to  the  strain  associated  with  an  applied  stress. 
Experimental  evidence  [94]  has  shown  that  the  mechanical  behavior  of  granular  media 
is  greatly  affected  by  their  anisotropy  which  is  closely  related  to  the  spatiad  arrange¬ 
ment  of  its  particles  and  the  fabric.  Knowing  the  mechanism  for  change  of  fabric 
during  deformation  will  provide  better  insight  into  the  evolving  anisotropy  of  gran¬ 
ular  materials.  Therefore,  the  general  concept  and  the  several  common  measures  of 
fabric  will  be  reviewed  in  the  following. 

The  precise  definition  of  a  fabric  tensor  varies  with  the  type  of  material 
and,  sometimes  for  the  same  material,  according  to  investigator.  The  choice  of  a 
particular  fabric  measure  is  a  matter  of  convenience  and  its  suitability  is  judged  by 
comparison  with  experimental  observation  [124].  A  relatively  universal  second-order 
moment  tensor  defined  by 

Nij  =<  riiTij  >  (2.1) 

where  <  ...  >  designates  the  sample  mean,  i.e.  <  n.-n^  >=  is  called  the 

anisotropy  tensor  by  Satake  (1982)  or  the  fabric  tensor  of  the  first  kind  by  Kanatani 
(1984).  In  eq.  2.1,  the  n,-  are  direction  cosines  of  the  cth  contact  normal  n  =  (n,)  to 
the  tangent  plane,  with  respect  to  the  orthogonal  coordinate  system.  Cn  is  the  total 
number  of  contax:ts  in  a  given  volume,  and  N  =  (iV,j)  is  symmetric  with  unit  trace. 

For  non-spherical  granules,  Nemat-Nasser  et  al.  (1983)  proposed  the  tensors: 

Hij  =<  m,m,  > 


(2.2) 


Hij  =<  mitij  >  (1.6) 

where  m,  is  the  Cartesian  components  of  a  unit  branch  vector,  a  branch  being  de¬ 
fined  as  the  connection  from  the  centroid  of  one  particle  to  that  of  another  touching 
particle.  They  even  suggested  the  inclusion  of  average  branch  length  /  and  contact 
area  a  into  the  fabric  tensor,  represented  by  E)q.  2.2  and  2.3,  to  reflect  additional  in¬ 
formation  on  the  microstructure.  Higher  order  fabric  tensors,  such  as  <  n.njTiin/  >, 
<  mifnjmkmi  >  and  <  riinjmkmi  >,  may  also  be  considered.  The  higher  order 
tensors  provides  more  information  regarding  the  details  of  the  anisotropy  (Kanatani 
1984  and  Subhash  et  al.  1991). 

Kanatani  (1984)  proposed  a  distribution  density  function  B(n),  defined  as 

£;(n)  =  T)Fij,...,kninj  . . .  n*  (2.4) 

in  which  t)  equals  to  1  /27r  for  two  dimensional  case  and  1  /Air  for  three  dimensional 
case,  a  tensor  of  even  rank  r,  is  referred  to  as  the  “rank  r  tensor  of  the  second 

kind”.  Then, 

£:(n)  =  £:(-n)  (2.5) 


/  Ein)dil  = 
Ja 


where  fl  is  the  surface  of  unit  sphere,  and  E{n)dCl  is  the  relative  number  of  normals  n 
falling  in  the  solid  angle  dU,  about  the  direction  n.  To  represent  the  density  function 
£(pl  by  the  second  rank  fabric  tensor  of  the  second  kind,  F,j,  we  have,  for  the  three 
dimensicnal  medium, 

£(n)  =  (2.7) 

Fii  =  -  i«iy)  (2.8) 

i,j  =  1,2,3,  and  for  a  two  dimensional  medium, 

^(n)  =  (2-9) 

Fii = mi  -  (2.10) 

for  t,j  =  1,2,  in  which  Su  denotes  Kronecker’s  delta. 
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2.2  Contact  Mechanics  and  Nonlinear  Elasticity 

Cohesionless  granular  materials  support  an  ambient  shear  stress  only  through 
the  contact  between  particles.  Therefore,  it  is  plausible  that  the  mechanism  of  local 
contacts  should  have  great  influence  on  overall  mechanical  properties  of  these  media. 
One  example  is  the  apparent  nonlinear  elasticity  at  small  strains  exhibited  collectively 
by  an  assemblage  of  particles  which  behave  individually  in  a  linear  elastic  way.  This 
effect  can  be  ascribed  to  the  intrinsic  nonlinearity  of  the  contact  mechanics  governing 
particle-particle  interactions  (Goddard  1990). 

Hertz  first  initiated  the  mathematical  study  of  the  effects  produced  by  mu¬ 
tual  compression  of  elastic  bodies  for  the  case  in  which  the  forces  between  bodies 
are  normal  to  the  contact  surfaces  [84,72].  Considering  two  elastic  spheres  in  contact, 
according  to  Hertzian  theory,  a  circular  contact  surface  is  produced,  with  radius  given 
by 

a  =  (2.11) 

where  /„  is  the  normal  force,  R  is  the  radius  of  the  spheres,  and  Mi  =  —  in 

which  E  and  v  are  Young’s  modulus  and  Poisson’s  ratio  of  the  material,  respectively. 
The  theory  also  gives  the  relative  approach  of  the  spheres 

6  =  2{MifjR}f^yf^  (2.12) 

Hence  the  apparent  normal  contact  stiffness  is  given  as 

The  tangential  stiffness  for  frictional  contacts  under  oblique  contact  force 
was  given  by  Mindlin  and  Deresiewicz  [84,26]: 

K  =  (2.14) 

where  M2  =  2(1  —  v)/{2  —  i/),  y  the  interpartide  friction  coefficient  and  ft  is  the 
resultant  shear  force  at  the  contact. 
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Figure  2.1:  Elastic  wave  velocity  in  an  FCC  packing  of  5  inch  diameter  steel  balls  with 
‘low’  (A)  and  ‘high’  (0)  dimensional  tolerances,  ±50  x  10"®  inches  and  ±10  x  10"* 
inches,  respectively;  after  Duffy  and  Mindlin  (1957).  The  broken  lines  of  slope  \ 
have  been  added  by  Goddard  (1990c)  for  comparison.  The  solid  lines  with  slope  5 
represent  the  Hertz-Mindlin  contact,  (a)  with,  and  (b)  without  tangential  stiffness. 
With  permission  of  the  author  of  [57]. 

The  above  theory  has  been  adopted  in  most  theoretical  treatments  of  the 
micromechanics  of  granular  media,  dating  from  the  landmark  w'orks  of  Mindlin  and 
co-workers  ([46]  etc.)  up  to  the  most  recent  publications  on  the  subject  ([44,134] 
and  the  references  cited).  As  suggested  by  above  equations,  the  underlying  theory 
leads  inevitably  to  the  pow-er-law  scaling  E  for  the  dependence  of  various  elastic 
moduli  on  confining  stress  p  and,  hence,  to  the  scaling  u  ~  pi  lor  various  elastic  wave 
velocities  v  (with  magnitudes  characterized  by  \jEJp,  here  p  is  the  material  density). 
However,  experimental  evidence  within  soil  mechanics  and  geophysics  show’s  that 
the  scaling  ~  pi  and  u  ~  pi  are  much  more  representative  (see  Goddard  [57] 
for  a  complete  survey),  although  the  pressure  dependence  may  change  from  p^^^  to 
p^l^  in  high  pressure  regime  or  under  prolonged  vibration  at  large  amplitudes.  Such 
observations  are  illustrated  in  figure  2.1. 

It  is  abo  reflected  quantitatively  in  the  widely  used  empirical  formula  for  the 
shear  modulus  of  dry  sands  (see  [57]  and  references  cited)  under  isotropic  confinement 
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at  initial  or  base  pressure  p: 


G  =  £l(«.  -  «)V{1  +  «))?'"  (215) 

where  e  is  void  ratio,  while  (  and  are  constants. 

In  a  detailed  analysis  [57],  Goddard  showed  that  one  can  explain  such  fre¬ 
quently  observed  departures  from  the  1 /6-power  dependence  predicted  by  Hertz- 
contact  theory  on  the  basis  of  two  rather  distinct  hypothesis.  The  first  involves 
nonhertzian  asperities  while  the  second  appeals  to  the  nonlinearities  arising  from 
strain-induced  changes  in  the  number  of  particle  contacts.  For  isotropic  confinement 
both  the  above  hypothesis  yield  a  1 /4-power  dependence  of  wave  speed  on  pressure 
at  low  confining  pressures,  with  a  transition  to  a  1 /6-power  dependence  at  high  pres¬ 
sures. 

2.3  Conduction  through  Granular  Materials 

2.3.1  Introduction 

The  prediction  of  the  effective  conductivity  of  two-phase  media,  in  which 
one  phase  is  dispersed  in  a  second,  h2is  occupied  engineers  and  physicists  for  the 
paist  one  hundred  and  twenty  years  [130,35,16,9,47,22,13,55].  This  long  interest  has 
been  fueled  by  the  proliferation  of  man-made  composite  materials  and  the  need  to 
predict  bulk  properties  such  as  effective  conductivity.  Maxwell  (1873)  was  the  first 
to  theoretically  calculate  the  effective  conductivity  of  a  dilute  stationary  suspension 
of  spherical  particles.  By  considering  only  the  interaction  of  a  single  sphere  in  a 
potential  gradient.  Maxwell  was  able  to  obtain  the  following  well-known  equation 

fc*  1  -f  2^<f> 

ko  I  — 

where,  ^  =  {a  —  !)/(« -1-  2);  a  is  the  ratio  of  conductivity  of  the  solid  particle  to  that 
of  the  matrix(or  fluid  phase);  k*  is  the  effective  conductivity  of  the  suspension;  ko  is 
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the  conductivity  of  the  matrix;  4>  is  solid  volume  fraction.  To  the  order  of  terms  in  (j) 
to  which  it  is  exact,  Eq.  2.16  takes  the  form: 

^  =  1  +  3)?^  +  0OT  (2.17) 

Kq 

A  hundred  years  later,  Jeffrey  (1973)  extended  Maxwell’s  result  to  by  the 

addition  of  two-sphere  interactions  for  a  random  haxd-sphere  dispersion 

T“  =  1  +  3/3^  -J-  0<f>^  -f  0{<f>^)  (2.18) 

Kq 

A  A 

where  —  ^{0)  is  a  slowly- convergent  infinite  series  in 

Some  progress  has  also  been  made  for  densely  packed  suspensions  of  perfectly 
conducting  spheres.  Keller  (1963)  solved  this  problem  correct  to  Inc,  where  e  is 
the  dimensionless  gap  width,  for  a  densely  packed  simple  cubic  array  of  spheres. 
Batchelor  and  O’Brien  (1977)  extended  Keller’s  work  to  include  touching  spheres 
and  near-perfect  conductors  by  using  a  mean-field  approach.  They  theoretically  treat 
the  thermal  or  electrical  conduction  through  static  particulate  media  in  the  limit  of 
maximum  volume  fraction,  for  which  the  particles  make  point  contact  with  each  other 
and  even  interface  with  flat,  convex  or  concave  surfaces  under  external  load.  They 
find  that  when  a  >>  1  the  effective  conductivity  of  random  two  phase  media  is  given 
by 

k” 

—  =4lna-ll  (2.19) 

Kq 

with  constant  4  predicted  by  the  theory  and  the  additive  constant  chosen  to  achieve 
a  reasonable  fit  with  a  variety  of  experimental  data  points.  Their  theory  suggests 
that  the  exact  method  of  forming  a  dense  suspension  will  strongly  affect  its  effective 
conductivity  by  the  resultant  average  coordination  niimber  of  the  particles,  and  it 
illustrates  that  the  microstructure  has  a  measurable  effect  on  the  conductivity  of  a 
suspension. 


2.3.2  Electric  Contacts 

In  present  study,  we  will  consider  a  simplified  two  phase  medium  with  the 
continuous  phase  being  nonconductive.  Electrical  conduction  through  a  packed  bed 
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of  steel  balls,  with  air  filling  the  interstitial  void,  represents  a  prime  example.  Under 
compressional  loading,  the  particles  are  pressed  together  and,  if  elastic,  will  deform 
slightly  and  will  develop  a  flat  circle  of  contact.  According  to  the  Hertz  theory 
described  in  section  2.2,  two  touching  elastic  particles  which  are  spherical  locally 
with  rsidius  R  will  develop  a  flat  contact  circle  of  radius 


(2.20) 


where  a  compression  force  /„  acts  on  each  particle  normal  to  the  common  tangent 
plane  at  the  point  of  contact.  Since  the  pore  fluid  is  non-conductive,  electric  flux  is 
only  possible  through  the  contact  circle,  and  the  distribution  of  potential  inside  the 
two  particles  is  approximately  the  same  as  that  of  the  velocity  potential  in  irrotational 
flow  of  an  incompressible  fluid  through  a  circular  hole  in  a  plane  wall  [16].  The  solution 
to  this  latter  problem  is  known,  and  shows  that  the  normal  flux  density  at  the  context 


circle  is 


j 

7r(a*  - 


(r<a) 


(2.21) 


and  the  total  current  across  the  circle  of  contact  is 


Q  =  f  J2wrdr  =  2akpA^  (2.22) 

Jo 

where  kp  is  the  conductivity  of  the  particle  material  and  A$  is  the  difference  in 
electrical  potential  between  the  particles.  From  the  above  equation,  the  contact 
resistance  between  two  particles  across  the  the  contact  circle  is 

2akp 

which  is  nothing  but  the  constriction  resistance  of  the  small  flat  contact  area.  A 
similar  equation  was  also  obtained  by  Holm  (1967)  and  Yovanovich  (1967). 

However,  the  contact  resistance  between  two  real  surface  is  far  more  compli¬ 
cated.  To  understand  why,  it  is  necessary  to  consider  the  nature  of  solid  surfaces  and 
the  dfect  of  foreign  materials  on  the  overall  resistance.  Contact  surfaces  are  irregular 
on  a  microscopic  scale.  Even  nominally  plane  surfaces  have  a  waviness  with  peak- 
to- valley  dimensions  typically  from  tenths  to  several  micrometers  [7,60,6].  When  two 
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contacts  are  brought  together  under  low  loads,  they  touch  at  only  a  few  asperities 
(multispot  contact).  As  the  load  is  increased,  more  asperities  come  into  contact  and 
the  surfaces  move  together.  Therefore,  the  true  area  of  contact  depends  on  normal 
load  and  the  hardness  of  the  materiad  [64,61].  This  zirea  is  only  a  small  fraction  of 
the  apparent  contact,  except  at  very  high  loads  where  the  surface  can  be  severely 
deformed.  Furthermore,  when  the  metal  surface  is  exposed  to  the  environment,  a 
contaminant  film  will  be  developed  through  processes  such  as  oxidation  and  corro¬ 
sion,  particulate  contamination  (airborne  and  wear  debris),  fretting,  etc.,  and  soon 
covers  up  the  virgin  metal  [136,6].  This  contaminant  film  is  often  extremely  noncon- 
ductive,  therefore,  preventing  electrical  conduction  through  the  contact.  Under  such 
circumstances,  conduction  is  not  possible  if  the  film  is  unbroken,  except  when  the 
film  is  only  a  few  atomic  layers  thick,  such  that  some  electron  current  can  penetrate 
it  by  means  of  the  tunnel  effect  [64].  Under  a  normal  mechanical  load,  the  insulating 
film  on  the  contact  asperities  deforms  plastically  and  fractures,  so  that  pure  metal 
substrates  are  once  again  exposed  to  each  other.  Figure  2.2  schematically  illustrates 
this  situation  where  the  apparent  area  of  contact,  the  metallic  regions,  and  the  places 
with  insulating  layers  are  differentiated.  The  lines  of  current  flow  converge  at  the 
region  of  metallic  contact,  called  “a”  spot,  as  illustrated  schematically  in  Figure  2.3. 
Contact  resistance  decreases  with  increasing  load.  The  softer  and  more  conductive 
the  metal,  the  lower  the  contact  resistance  will  be  at  a  given  force. 

2.3.3  Multispot  Theory  of  Contact 

Generally,  a  multispot  problem  is  simplified  by  assuming  all  of  the  a- spots 
to  be  circular  and  to  lie  at  distances  from  each  other  which  are  large  compared  to  the 
radii,  thus  permitting  the  assumption  of  no  interference  between  different  a-spots. 
Thus,  the  total  resistance  becomes 

^  2k„E(ii 

where  subscript  ‘t’  represents  the  I'th  a-spot. 
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Figure  2.2:  Schematic  illustration  of  apparent  contact  surface.  The  metallic  contact 
regions  a  are  indicated  by  dark  areas.  Contact  at  region  h  (shaded  areas)  is  with 
insulating  contaminant  film.  Region  c  does  not  touch. 


solid  surfaces 


Figure  2.3:  Microscopic  view  of  a  real  contact  interface:  “a”  spot  and  lines  of  current 
flow. 
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For  the  case  in  which  the  a-spots  lie  close  to  each  other  so  that  the  constricted 
lines  of  flow  from  different  a-spots  deflect  each  other,  then  Eq.  2.24  is  no  longer  valid. 
Holm  [64]  has  made  some  approximations  for  the  case  of  the  uniformly  distributed 
a-spots,  giving  the  following  expression 

1 


i2c  = 


2irnak, 


■arctan— - 0.6-^^ — : - h 


kpA-r 


ikpr 


(2.25) 


where  n  is  number  of  a-spots,  a  the  radius  of  a-spot,  21  the  average  distance  between 
neighboring  a-spots,  Ar  the  area  of  apparent  contact,  and  r  the  radius  of  the  apparent 
contact  surface. 


2.3.4  Effective  Conductivity 

Suppose  that  a  uniform  intensity  gradient  is  set  up  in  the  medium,  perhaps 
by  imposing  uniform  and  different  values  of  the  intensity  at  two  distant  parallel 
boundaries.  Although  we  restrict  ourselves  to  the  electrical  conduction  problem, 
the  formulation  can  be  applied  to  the  transport  of  other  scalar  properties  such  as 
thermal  conduction  and  mass  diffusion.  Henceforth  we  shall  use  terms  and  notation 
appropriate  to  the  case  of  electrical  conduction  for  convenience.  So  the  mean  intensity 
gradient  will  be  written  as  <  >,  where  is  the  electrical  potential  gradient  at 

a  point  in  the  medium  (not  necessarily  lying  in  the  matrix)  and  the  brackets  <  ...  > 
denote  an  average  over  the  entire  volume  of  the  medium.  The  local  current  density 
J  is  equal  to  —koV^  at  a  point  in  the  matrix  and  —kpV^  at  a  point  in  a  particle. 
At  each  point  on  the  surface  of  a  particle  $  and  the  normal  component  of  J  are 
continuous;  and  at  each  point  not  on  such  a  surface 

V  •  J  =  0  and  =  0.  (2.26) 

Because  of  the  intrinsic  linearity,  the  magnitude  of  all  potential  differences 
are  proportional  to  the  magnitude  of  <  >,  and  so  for  the  mean  flux  density  we 

have  the  linear  relation  [16] 


<  J  >=  -K*  <  >, 


(2.27) 
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where  the  effective  conductivity  K*  is  a  second- rank  tensor,  dependent  on  the  struc¬ 
ture  of  the  medium. 

Next,  we  will  derive  a  formulation  for  the  mean  flux  density  through  gran¬ 
ular  media  with  a  nonconductive  fluid  phase.  In  this  case,  if  particles  are  highly 
conductive,  the  resistivity  will  be  brought  about  at  contact  points  only.  Therefore, 
the  medium  can  be  approximated  by  a  resistor  network,  as  shown  in  figure  2.4,  with 
particle  centers  being  nodes  and  interparticle  contact  being  resistor.  Similarly,  Fig.  2.4 
can  also  be  used  to  represent  a  elastic  network  [57]  by  replacing  resistors  by  elastic 
springs  whenever  elastic  properties  are  involved.  By  definition,  we  have 

(2.28) 

or 

(2.29) 

where  V,  Vp  and  Ki  represent  the  total  volumes  of  medium,  particle  and  fluid  phase, 
respectively,  N  the  number  of  particles. 

For  the  nonconductive  matrix,  the  first  term  on  the  right  in  Eq.  2.29  vanishes, 
so  we  have 

(2.30) 

Applying  Gauss’  divergence  theorem  to  Eq.  2.30,  the  following  is  obtained 

1  ^  r 

<  J  >=  -  xJ  •  ndS  (2.31) 

where  Sp  is  the  surface  of  a  particle,  x  is  the  spatial  position  of  the  points  on  Sp,  and 
n  represents  the  unit  normal  to  Sp.  Noting  that  the  term  J  •  ndS  is  current  passing 
through  the  portion  dS  of  particle  surface  Sp  and  assuming  that  the  current  passes 
through  the  individual  particle  at  discrete  points  of  contact,  we  obtain 

<j>=if;f;x(3  (2.32) 

where  Q  is  current  flowing  through  the  points  of  contact.  For  spheric2d  particles, 
X  =  x**  -j-  /2n,  here  x**  is  the  position  vector  of  the  particle  centroid,  and  Eq.  2.32 


r 


I 


I 

I 


I 


17 


Figure  2.4:  Network 


becomes 

where,  electrical  conservation  (Kirchhoff’s  law)  gives: 


5:x'’g=x»>5:Q=o 


(2.33) 


(2.34) 


The  formulation  described  above  is  similar  to  one  to  derive  the  macroscopic  stress 
tensor  for  a  assemblage  of  granular  materials  except  for  the  tensorial  orders  involved. 


In  order  to  compute  the  interparticle  current  Q,  the  local  potentials  $  of  the 
particles  must  be  determined  first.  These  local  potentials  can  be  further  divided  into 
two  parts:  one  derived  from  the  mean  potential  gradient,  the  other  being  fluctuation 
necessary  to  satisfy  the  current  balance  condition  within  the  system.  The  latter  is 
obtained  by  solving  the  system  of  linear  equations 

A$'  =  B  (2.35) 

where  A  represents  the  conductance  matrix,  is  the  fluctuation  and  B  is  the  net  un¬ 
balanced  current  owing  to  the  mean  potential  gradient.  When  a  cluster  of  one  or  more 
particles  is  isolated  from  the  rest,  the  matrix  A  become  singular.  This  singularity  is 
resolved  by  means  of  the  relaxation  method  to  be  discussed  in  Section  3.1.3. 


From  Eq.  2.27  and  Eq.  2.34,  one  can  infer  the  effective  conductivity  tensor 
K*,  which  represents  exact  solution  of  the  problem,  in  contrast  to  the  mean-field 
theory  of  Batchelor  &  O’Brien  [16]  described  next. 

Their  mean- field  theory  assumes  that  the  potential  difference  $j— between 
particles  i  and  j  is  given  by  (x^  -  x*)-  <  >=  — 2i2n*  <  >,  i.e.  to  the 

difference  between  the  potential  at  the  two  sphere  centers  in  the  mean  potential 
field.  Furthermore,  we  assume  that  the  contact  resistances  at  all  contacts  take  the 
identical  average  value  Rc.  With  these  assumptions,  one  no  longer  has  local  electrical 
conservation,  but  rather  a  global  conservation  in  some  average  sense. 

According  to  above  assumptions,  we  can  write 


_  1  .  .  2R 

Q  =  -g-A$  =  — g-n-  <  > 

Rc  Rc 


(2.36) 


Combining  Eq.  2.36  with  Eq.  2.33  leads  to  the  result  similar  to  one  given  by  Batchelor 
&  O’Brien 


<  J  >= 


2jR2  ® 

<  V$> 


R^V 


(2.37) 


Therefore,  by  comparing  Eq.  2.37  with  2.27,  one  obtains  the  effective  conductivity 
tensor  according  to  the  mean-field  theory 


2/32  N  c 


RcV 


(2.38) 


or  in  terms  of  the  fabric  tensor  N, 


2R?Cn 


(2.39) 


where  Cjf  is  the  total  number  of  contacts  in  a  given  volume,  which  is  equal  to  N 
times  the  average  coordination  number. 

The  estimates  of  effective  conductivity  provided  by  the  mean-field  theory 
of  Batchelor  &  O’Brien,  as  outlined  in  [16],  is  different  from  that  provided  by  the 
conventional  Voigt- Reuss-Hill  bounds  which  do  not  depend  on  the  fabric  tensor. 


Chapter  3 


Quasi-static  Simulation  and 
Assemblage  Generation 


As  a  continuation  and  extension  to  three  dimensions  of  the  work  of  Bashir 
and  Goddard  [14],  we  have  developed  an  improved  version  of  their  2D  programs 
by  introducing  several  new  techniques  [58].  Among  these  is  the  combination  of  the 
particle-assemblage  generation  and  the  computation  of  particle  motion  together  into 
one  single  program.  The  new  program  includes  a  shuffling  algorithm,  for  generating  an 
initially  random  loose-packed  configuration  of  particles,  and  an  improved  microcell- 
adjacency  method  to  further  accelerate  particle-contact  search.  Furthermore,  we  have 
also  overcome  a  singularity  in  the  quasi-linear  system  of  equilibrium  equations  by 
means  of  a  relaxation  method[121].  Our  program  is  able  to  simulate  any  deformation 
history  and  allows  us  to  study  both  mechanical  and  scalar  transport  properties  of  an 
idealized  granular  assemblage  simultaneously. 

The  validity  of  the  numerical  algorithm  is  “tested”  by  comparison  against 
the  triaxial  compression  experiments.  I  know  of  no  other  way  to  test  it  except  against 
other  numerical  codes,  however,  which  have  their  own  problems.  The  triaxial  com¬ 
pression  experiments  are  to  be  described  in  Chapter  4,  while  the  qualitative  compar¬ 
isons  between  the  results  of  numerical  simulations  and  triaxial  compression  tests  will 
be  given  in  Section  6.2. 


20 


The  following  sections  will  describe  the  quasi-static  simulation  and  afore¬ 
mentioned  techniques. 

3.1  The  Model  and  The  Relaxation  Method 

Unlike  dynamic  simulations,  in  which  the  full  (Newtonian)  dynamical  equa¬ 
tions  are  employed  to  update  particle  configurations,  the  present  quasi-static  scheme 
first  moves  every  particle  in  the  system  according  to  the  mean  deformation  gradient, 
thus  destroying  the  state  of  equilibrium.  Hence,  particles  have  to  be  relocated  to 
a  equilibrium  position  by  means  of  fluctuations  about  the  mean.  These  fluctuating 
displacements  of  an  individual  particle  are  determined  by  the  total  unbalanced  elas¬ 
tic  force  exerted  on  it  as  a  result  of  the  mean  deformation.  Equilibrium  is  achieved 
by  an  algorithm  that  allows  the  system  to  expand  or  contract  volumetrically  when 
necessary  to  maintain  a  control  pressure  or  stress  at  a  desired  level  which,  thereby, 
allows  us  to  compute  the  granular  dilatancy. 

3.1.1  The  Force-displacement  Law 

During  the  deformation  of  granular  assemblages,  particles  move  with  inde¬ 
pendent  degrees  of  freedom  and  interact  with  each  other  only  at  their  contact  points. 
The  assumed  force-displacement  relationship  will  be  presented  here  for  the  case  of 
two  spherical  particles  A  and  B  in  contact,  as  shown  in  figure  3.1. 

Particle  radius  is  denoted  by  R  and  its  centroid  by  X.  Upon  the  deformation, 
a  particle  undergoes  translational  and  rotational  displacement  increments  u  and  w, 
respectively.  The  superscript  in  Fig.  3.1  and  in  the  sequel  denotes  a  given  particle. 
The  unit  contact  normal  vector  to  the  tangential  plane,  viewed  from  A  to  B,  is 
expressed  as  n  =  (X^  —  X^)/|X^  —  X^|.  The  interaction  between  the  particles 
depends  on  the  relative  motion  of  the  contact  points.  The  vectorial  components  of 
relative  displacement  in  normal  and  tangential  direction  are  written  as: 


Aun  =  (u^  -  u"*)  •  nn 


(3.1) 


Figure  3.1:  Interaction  between  two  particles 


and 

Aut  =  (u^  -  -  Au«  +  R®  X  X  (3.2) 

where  ==  u®  -  u"*,  R^  =  and  R®  =  R^n. 

These  relative  displacements  are  used  to  calculate  increments  of  normal  and 
shear  forces,  Afn  and  Aft,  according  to: 

Afn  =  fc„Aun  (3.3) 

and 

Af,  =  ktAut  (3.4) 

where  kn  and  kt  denote  the  normal  and  tangential  elastic  stiffnesses,  respectively, 

which  may  be  allowed  to  depend  on  Afn  and  Aft.  However,  since  we  are  primarily 

interested  in  nearly  rigid  particles,  the  exact  dependence  of  the  elastic  stiffnesses  on 
Af  is  presumably  not  important  (vide  infra). 

Furthermore,  the  force  increments  Afn  and  Aft  are  added,  respectively,  to 
the  forces  f®  and  f?  that  existed  previously  between  two  particles  to  yield  the  current 
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values: 

f„  =  f^  +  Af„ 

and 

f,  =  Jf  +  Af, 


(3.5) 

(3.6) 


The  components,  f„  =  /„n  and  f«  =  /tt,  of  the  force  vector  act  along 
the  directions  of  contact  normal  and  tangent  plane  and  are  both  set  to  be  zero  if 
/n  is  not  compressions!  (since  cohesionless  particles  cannot  sustain  a  tensile  force). 
A  (Coulomb)  sliding  friction  law  is  incorporated  as  follows:  the  magnitude  of  the 
shear  force  ft  given  by  eq.  3.6  is  checked  against  the  maximum  possible  shear  force 
magnitude: 

(/t)mo®  =  |fn|  "I"  Oh  C^'^) 

where  fi{=  tan  <f>^)  is  the  coefficient  of  sliding  friction  (defining  the  so-called  2aigle  of 
intergranular  friction  ^^),  and  Ch  represents  cohesion,  which  is  taken  as  identically 
zero  for  the  non-cohesive  particles  considered  here.  If  |ft|  exceeds  (/t)mox>  sliding 
occurs  at  the  contact  point.  Under  this  circumstance,  ft  takes  the  value  of  (/t)mo*» 
and  maintains  its  direction.  Therefore,  the  total  force  and  couple  exerted  on  particle 
A  by  particle  B  are  given  by: 


f  =  fn  +  ft 

(3.8) 

M  =  X  f, 

(3.9) 

3.1.2  The  Governing  Equations 

The  force  f  and  couple  M  are  next  decomposed  into  three  Cartesian  com¬ 
ponents,  which  yields  in  matrix  form: 

F  =  F^  +  k^BAu(AB)  (3.10) 

where,  F  =  [/*,  /v,  /n  wirj  ”**]»  generalized  force,  represents  the  components  of 
force  and  moment  exerted  currently  on  particle  A  by  B.  F°  =  [/°,  ,  m°,  m°,  m®] 

represents  the  components  of  the  force  and  the  moment  in  the  previous  state.  The 


matrix  kxB  is  called  the  local  contact-stiffness  matrix,  while  Au{AB)  is  the  general¬ 
ized  relative  displacement  between  A  and  B,  written  as: 


Au{AB)  = 


uf-ui 
uf-uf 
uf  -  uf 


-I- 

r^oJy  + 


\  -I-  j 


(3.11) 


where  the  subscripts  denote  the  corresponding  Cartesian  component. 

For  all  contacts  on  A  to  be  in  static  equilibrium,  the  sum  of  generalized  force 
must  vanish: 

5;F  =  5;F»+5;k^BAu(AB)  =  o  (3.12) 

B  B  B 

or 

£kABAu(AB)=-Z^  (3.13) 

B  B 

In  the  current  simulation,  the  displacement  of  each  particle  is  additively 
decomposed  into  two  components:  the  macroscopically  imposed  mean  IT  defined  by 
the  global  velocity  gradient  and  a  fluctuation  u',  the  latter  being  such  that  the  force 
balance  eq.  3.13  is  satisfied.  Therefore,  eq.  3.13  becomes: 


Y,  Au'(AB)  =  -  53  F°  -  5;  k„B ATr( AB)  (3.14) 

B  B  B 

for  A  =  1, 2,  •  • . ,  iV,  with  N  denoting  the  total  number  of  particles  within  the  system. 
This  represents  a  system  of  quasi-linear  equations  for  the  {6N)  fluctuating  particle 
displacements: 


Kx  =  b 


(3.15) 


where  K  is  the  grand  stiffness  matrix,  x  =  {u^(l),Uy(l),u'(l),u;i(l),u>'(l),u;'(l),  •  •  *, 
ui(iV),  Uy{N),  u' (JV),  w' (JV),  w' (iV)]  the  vector  of  the  fluctuating  displacements 

and  rotations,  and  b  the  unbalanced  force  arising  from  the  mean  displacements  and 
forces  from  the  prior  deformation  step. 
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3.1.3  Relaxation  Method 

The  K  matrix  in  (  3.15)  becomes  singular  whenever  a  cluster  of  particles 
is  isolated  from  the  remainder,  giving  rise  to  “neutral”  or  “zero-frequency”  elastic 
vibrational  modes  representing  a  finite  dimensional  null  space  of  K.  To  overcome 
this  singularity,  Bashir  and  Goddard  add  a  small  artificial  term  to  the  diagonal  of  K 
and  then  employ  a  gaussian  elimination  method  to  solve  (3.15)  for  x. 

In  the  present  simulation,  we  utilize  the  relaxation  method  (origin<dly  due 
to  Southwell  [121])  as  our  linear-equation  solver,  since  it  effectively  cuts  out  the  zero- 
frequency  modes  of  K.  Being  an  iterative  method,  relaxation  involves  two  procedures 
to  accelerate  convergence.  First,  the  relaxation  order  is  determined  by  searching 
for  the  residual  of  greatest  magnitude  |i2,-|max  (the  residuals  R  being  the  difference 
between  the  right-hand  and  left-hand  sides  of  (3.15),  evaluated  at  the  current  values 
of  X  in  an  iteration),  then  “relaxing”  the  corresponding  equation  by  calculating  a 
new  value  of  Xi  so  that  {Ri)mas  —  0.  This  modifies  all  other  residuals,  which  also 
depend  on  Xj.  The  procedure  is  applied  repetitively  until  all  the  residuals  satisfy  a 
preset  convergence  criterion  on  some  norm  |f2|.  In  the  present  context  the  fluctuations 
determined  by  the  relaxation  method  serve  to  move  only  those  particles,  or  particle 
clusters,  having  non-equilibrated  forces  or  moments.  Hence,  isolated  clusters  do  not 
fluctuate,  and  we  avoid  the  singularity  in  inverting  (3.15). 

For  the  packing  algorithm  described  in  the  sequel,  the  relaxation  method 
is  particularly  effective,  since  in  the  early  stages,  the  number  of  particle  contacts  is 
small  and  only  those  particles  not  in  equilibrium  need  be  moved.  Furthermore,  the 
relaxation  scheme  always  finds  the  maximum  unbalanced  forces  and  adjusts  particle 
positions  so  as  to  balance  them  out. 
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1  2  ...  nm 


(a) 


(b) 


Figure  3.2:  Microcells  and  the  simulation  cell  in  (a)  initial,  and  (b)  sheared  configu¬ 
rations 

3.2  The  Microcell  Method  and  The  Adjacency 
Matrix 

In  the  computer  simulation  of  a  classical  mechanical  system  of  N  interacting 
particles,  it  is  generally  necessary  to  search  for  all  particles  within  the  range  of  spatial 
interaction  of  a  given  particle.  In  general,  one  needs  N{N  —  l)/2  such  searches, 
including  a  time-consuming  evaluation  of  particle  separations,  a  non- trivial  task  when 
the  number  of  particles  is  large.  However,  the  search  time  can  be  reduced  to  0{N)  by 
means  of  spatial  microcell  methods  [3]  and  the  associated  adjacency-matrix  technique. 

In  the  2D  case,  for  instance,  the  deformable  simulation  cell  is  divided  into 
regular  lattice  of  nm  x  m  initially  square  microcells  as  shown  in  Fig.  3.2.  A  microcell  is 
small  enough  to  contain  the  center  of  at  most  one  particle  throughout  the  subsequent 
cellular  deformations.  All  microcells  are  then  labded  ordinally.  For  each  microcell, 
the  definition  of  adjacent  microcells  may  include  a  neighborhood  extending  several 
microcell  layers  outward,  depending  upon  the  range  of  the  pair  interaction  considered. 
Whatever  the  range,  a  matrix  Ac  defines  the  adjacency  of  microcells: 


Ac{i,j)  = 


(3.16) 


1 ,  if  microcells  t  and  j  are  adjacent 

« 

0,  otherwise 

where,  i,j  =  1,2, ...,  mc{mc  =  nm  x  nm),  which  is  nothing  more  than  the  connectivity 
matrix  of  the  associated  graph  [131].  We  next  define  a  second  matrix  Oc  to  represent 
the  occupancy  of  microcells  by  particles,  such  that: 

il,  if  microcell  i  is  occupied  by  particle  j 

(3.17) 

0,  otherwise 

where,*  =  1,2, ...,  me,  and  j  =  1,2,  ...,iV.  A  third  matrix  j4p  is  then  used  to  represent 
the  adjacency  of  two  particles: 

{1 ,  if  particles  t  and  j  are  adjacent 

(3.18) 

0,  otherwise 

where,  i,j  =  1,2, ...,  N,  as  determined  by  their  occupancy  of  adjacent  or  non-adjacent 
microcells.  This  matrix  gives  the  “Verlet  neighbor  list”  of  molecular  dynamics  [3]  and 
can  be  expressed  as  the  matrix  product: 

Ap  =  AcOc  (3.19) 

Once  the  microcell  adjacency  matrix  Ac  is  established,  it  remains  unchanged 
as  long  as  the  microcell  topology  is  fixed  during  the  simulation.  Upon  determination  of 
the  occupancy  matrix  Oc  at  each  deformation  step,  the  particle  adjacency  matrix  Ap 
can  be  found  easily  by  the  simple  operation  (3.19).  However,  (3.19)  is  computationally 
equivalent  to: 

Ap{i,j)  =  Ac{map{i),map{j))  (3.20) 

where  map  defines  a  mapping  array  whose  element  map{i)  equals  the  ordinal  number 
(1,2, ...,  me)  of  the  microcell  occupied  by  particle  i  and  which,  therefore,  corresponds 
to  the  row  vectors  of  Oc.  Based  on  the  computed  particle  adjacency  matrix  Ap,  the 
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current  program  searches  for  all  particle  contacts  in  order  to  construct  the  stiffness 
matrix  and  to  calculate  contact  forces. 

In  our  simulation  on  rigid  spheres,  the  size  of  the  cubical  microcell  is  chosen 
such  that  its  largest  diagonal  is  equal  to  the  smallest  particle  diameter  in  the  system 
to  assure  that  not  more  than  one  particle  simultaneously  occupies  a  given  micro¬ 
cell.  Furthermore,  the  largest  particle  diameter  defines  a  cutoff  distance  at  which  one 
must  search  for  a  potential  contact  with  a  neighboring  particle.  For  a  2D  monodis- 
perse  disk  assemblage,  therefore,  two  surrounding  layers  will  be  sufficient  to  cover 
the  cutoff  distance,  which  means  there  are  24  microcells  adjacent  to  each  microcell. 
At  the  start  of  a  simulation,  the  microcell  adjacency  matrix  Ac  is  constructed  ac¬ 
cordingly  and  remains  unaltered  throughout  the  computation.  At  each  deformation 
step,  whenever  particles  move  to  new  positions  the  mapping  array  is  updated,  which 
is  a  rapid  process.  For  a  given  particle,  we  need  only  look  at  the  neighboring  24 
microcells  surrounding  its  microccii  to  find  the  adjacent  particles.  In  the  worst  case, 
24  searches  would  be  required  if  all  neighboring  microcells  were  occupied.  Therefore, 
12iV  provides  an  upper-bound  on  the  total  searches  necessary  if  we  consider  a  pair  of 
neighboring  particles  as  one  search. 

In  reality,  the  number  of  searches  required  depends  upon  the  number  of 
particles  lying  within  the  cutoff  distance  or  upon  the  system  density  and  configuration. 
In  a  2D  random  particle  assemblage,  the  average  number  of  necessary  searches  is  far 
less  than  24  per  particle.  Through  our  computations,  we  have  found  that  the  average 
number  of  searches  for  each  particle  is  about  6  for  random  loose-packed,  and  1 1  for 
random  dense-packed  monodisperse  disk  systems.  Therefore,  the  total  number  of 
searches  is  approximately  3N  and  5.5iV,  respectively. 


3.3  Random  Configuration  Generation 


In  the  past  thirty  years  or  so,  the  packing  of  disks  and  spheres  of  equal 
radii  in  2D  and  3D  has  been  studied  extensively  by  both  experimental  and  theoret- 
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ical  means,  in  part  because  they  serve  as  useful  models  for  a  variety  of  amorphous 
materi2ds  such  as  molecular  fluids  and  glasses.  The  macroscopic  properties  of  gran¬ 
ular  materials  and  porous  media  have  also  been  modeled  by  sphere  packings.  Three 
models  which  are  commonly  employed  for  packings  of  disks  and  spheres  are  the  dense 
ordered  packing,  dense  random  packing,  and  loose  random  packing.  The  dense  or¬ 
dered  packing  for  rigid  spheres  of  equal  radii  occurs  when  the  density  is  equal  to  0.7405 
in  3D(FCC  or  HCP).  Similarly,  the  density  is  equal  to  0.9069  in  2D(triangular).  For 
dense  random  packings,  it  is  generally  believed  that  the  densities  fall  into  a  range 
0.62  to  0.66  for  3D  and  0.81  to  0.87  for  2D  [19,50,71,114,133]. 

In  a  previous  study  of  2D  disk  assemblages,  Bashir  and  Goddard  [14]  have 
employed  two  distinct  algorithms  to  generate  two  types  of  assemblages:  imperfect  tri¬ 
angular  close-packed  for  the  monodisperse  assemblage  and  pseudogravitational  for  the 
polydisperse.  Recognizing  the  limitation  of  those  algorithms  in  allowing  for  variable 
initial  densities  and  for  generating  random  isotropic  configurations,  we  have  developed 
a  packing  algorithm  which  is  capable  of  densifying  an  initially  random  loose  config¬ 
uration  to  any  desired  density  for  both  monodisperse  and  polydisperse  assemblages 
by  means  of  cyclic  shear  under  isotropic  confining  pressure.  One  could  if  desired  add 
body  forces  such  as  gravity  but  we  shall  not  consider  them  here. 

3.3.1  The  Shuffling  Algorithm 

There  are  many  ways  of  realizing  random  sequences,  the  conventional  one 
being  the  standard  random- number  generation.  For  the  purpose  of  generating  random 
particle  assemblages,  we  introduce  a  new  way  of  rapidly  obtaining  repeated  random 
sequences  of  numbers  by  means  of  a  card-shulfling  algorithm.  The  idea  is  to  degrade 
the  order  of  a  given  set  of  numbers  by  means  of  a  certain  number  of  “riffle”  shuffles. 

In  shuffling  theory  [43],  a  single  riflBe  shuffle  can  generate  at  most  two  in¬ 
creasing  sequences  for  an  ordered  n-member  set  5,  where  a  increasing  sequence  is 
defined  to  be  a  sequence  whose  members  are  in  the  increzising  order  of  their  ordinal 
numbers  in  the  original  set  5.  If  F’n(i7)  be  the  number  of  permutations  of  n  items 
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with  exactly  R  increasing  sequences,  it  can  be  shown  that  [108]: 


f.(«)  =  D-iW  "  ‘ 

i=0  j 


(3.21) 


Shannon’s  theorem  [51]  states  that  a  permutation  with  exactly  2*  increasing  sequences 
can  be  obtained  by  k  riffle  shuffles  of  the  original  set  members  in  only  one  way. 

Hence,  in  k  riffle  shuffles,  the  total  number  of  permutations  that  can  be 
achieved  is:  ^ 

Uk)  =  i:  F4R)  (3.22) 


R=1 


and  thus,  the  number  of  riffle  shuffles  ib,  sufficient  to  generate  a  sequence  of  random 
numbers  is  given  by  [43]: 


T^{k)  >  n! 


(3.23) 


(3.23)  implies  that  for  a  deck  of  n  =  52  cards,  fc  =  7  riffle  shuffles  are  sufficient  to 
obtain  a  nearly  random  sequence  [43]. 

In  order  to  compare  our  shuffling  algorithm  with  a  standard  random-number 
generator,  we  have  calculated  the  autocorrelation  between  the  shuffled  sequence  of 
fifty  numbers  and  an  initial  ordered  sequence,  as  well  as  the  autocorrelation  between 
two  random  sequences  obtained  by  a  random  number  generator.  Let  S{i)  represent 
the  elements,  i  —  1,2,  ...,n,  of  an  n-member  sequence,  then,  we  employ  as  autocorre¬ 
lation  function  between  two  such  sequences  and  the  formula: 


which  treats  the  sequences  as  cyclical.  We  have  also  measured  the  CPU  time  required 
by  both  techniques,  and  the  details  will  be  discussed  below. 


3.3.2  Shuffling  vs  The  Random  Number  Generator 

In  our  shuffling  algorithm,  a  variant  of  the  riffle  shuffle  is  used,  wherein  each 
shuffle  consists  of  one  random  “cut”  and  “flip”,  and  one  interlacing  shuffle  [43].  To 
compare  our  shuffling  algorithm  with  a  random  number  generator,  we  have  computed 
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Figure  3.3:  A  comparison  between  the  shuffling  algorithm  and  random  number  gen¬ 
erator 

the  autocorrelations  between  the  shuffled  sequences  of  fifty  numbers  and  the  initial 
ordinal  sequence  as  a  function  of  number  of  shufflings.  Similarly  we  also  computed 
the  autocorrelations  between  the  first,  second  and  succeeding  sequences  with  the  first 
sequence  of  fifty  random  numbers  generated  by  a  random  number  generator  in  our 
computer(a  HP-730^^  workstation).  Fig.  3.3  indicates  that  sequences  of  random 
number  generated  by  the  shuffling  algorithm  are  as  good  as  those  obtained  by  the 
random  number  generator  in  terms  of  randomness,  even  when  we  employ  only  seven 
riffle  shuffles  for  all  of  our  computations.  To  compare  performance  in  speed,  we  have 
measured  the  CPU  time  required  by  both  techniques.  Fig.  3.4  clearly  shows  that  the 
random  number  generator  uses  approximately  four  times  as  much  CPU  time  as  the 
riffle  shuffle. 

3.3.3  Initial  Random  Loose  Configurations 

In  our  packing  algorithm  the  size  of  microcell  is  chosen  sufficiently  small  so 
as  to  contain  the  center  of  not  more  than  one  single  particle  under  any  circumstance. 
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Figure  3.4:  A  comparison  of  speeds  between  the  shuffling  algorithm  and  random 
number  generator 

but  also  sufflciently  large  so  as  to  minimize  the  total  number  of  microcells. 

Fig.  3.5  shows  a  2D  rectangular  microcell  ABCDA  with  two  sides  being 
denoted  by  Ax  and  Ay.  When  subject  to  a  simple  shear,  the  microcell  ABCDA  is 
deformed  to  ABC'D'A.  If  the  largest  dimension  AC  is  chosen  to  be  equal  to  the 
smallest  particle  diameter  we  are  assured  that  no  two  particles  can  simultaneously 
occupy  the  same  microcell  throughout  the  deformation.  Therefore,  we  have: 
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Figure  3.5:  Microcell  geometry 
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Ay  =  -  (3.25) 

Y  ('Ymo*  "I"  ^xy)  4"  1 

Ax  =  Ay  •  Vxy  (3.26) 

where,  dmin  is  the  smallest  particle  diameter,  7moe  is  the  magnitude  of  the  total  shear 
strain,  and  r^y  is  the  ratio  of  two  sides  (r^y  =  1,  usually). 

By  assuming  an  initial  density  and  knowing  the  total  particle  volume,  we 
estimate  the  size  of  the  simulation  cell  and  divide  the  cell  into  a  lattice  of  nm  x 
nm  microcells.  The  microcells  are  then  labeled  ordinally  from  1  to  me  (=nm  x 
nm).  To  place  N  particles  randomly  in  the  simulation  cell,  we  generate  a  random 
sequence  of  microcells  by  employing  the  shulBing  algorithm  described  before.  We 
then  pick  a  microcell  from  the  random  sequence  and  place  a  particle  randomly  within 
the  microcell,  so  as  to  avoid  overlap  with  previously  placed  particles,  until  all  N 
particles  are  placed  successfully.  Otherwise,  the  number  of  microcells  is  increased 
and  whole  process  is  repeated.  Fig.  3.6  and  3.7  show  one  such  random  loose-p^^cked 
configuration  for  132  disks  and  its  corresponding  random  dense-packed  configuration. 
Similar  configurations  for  48  polydisperse(multiple-sized)  spheres  are  displayed  in 
Fig.  3.8  and  3.9. 

3.3.4  Radial  Distribution  Functions 

The  radial  distribution  functions  g{r)  for  the  monodisperse  assemblages,  r 
being  scaled  on  particle  diameter,  have  been  computed  and  compared  with  those  gen¬ 
erated  by  the  Percus-Yevick  (P-Y)  equation  of  statistical  mechanics  and  by  a  Monte 
Carlo  (M-C)  simulation  [127,132,137],  to  verify  that  both  our  loose-packed  configura¬ 
tion  and  dense-packed  systems  are  random  for  2D  as  well  as  3D  assemblages.  Fig.  3.10 
shows  the  smoothed  g{r)  distribution  function  for  100  realizations  of  an  initially  loose- 
packed  configuration  of  132  disks(density=0.43).  Fig.  3.11  shows  the  smoothed  g(r) 
for  100  retdizations  of  dense-packed  configurations  of  132  disks(density=0.80,  close  to 
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Figure  3.9:  Random  dense-packed  configuration  and  associated  contact  bond  net¬ 
work  for  48  poly-disperse  spheres(density=0.60).  The  thickness  of  rods  represents 
the  scaled  magnitude  of  normal  force  between  particles. 
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Figure  3.10:  The  simulated  radial  distribution  function  for  2D  loose-packed  config¬ 
urations  of  132  disks(discrete  points),  compared  with  the  P-Y  radial  distribution 
function(solid  curve) 

those  for  2D  random  dense  packing).  For  3D,  Fig.  3.12  shows  the  smoothed  ^(r)  for 
30  realizations  of  loose-pax:ked  configurations  of  132  spheres(density=0.27).  Finally, 
Fig.  3.13  shows  g{r)  for  10  realizations  of  moderately  dense-packed  configurations 
of  90  spheres(density=0.58).  Our  computed  radial  distribution  functions  reveal  that 
both  the  initially  loose-packed  and  the  final  dense-packed  systems  are  quite  random, 
at  least  if  one  accepts  the  molecular  model  as  a  standard. 

The  text  of  this  chapter,  in  part,  appeared  in  [58].  The  dissertation  author 
was  the  secondary  author  of  the  publication,  and  shared  equal  responsibility  with 
other  co-authors. 
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Figure  3.11:  The  simulated  radial  distribution  function  for  2D  dense- packed  config¬ 
urations  of  132  disks(discrete  points),  compared  with  the  M-C  radial  distribution 
function(solid  curve) 


Figure  3.12:  The  simulated  radial  distribution  function  for  3D  loose-packed  config¬ 
urations  of  132  spheres(discrete  points),  compared  with  the  P-Y  radial  distribution 
function(solid  curve) 
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Figure  3.13:  The  simulated  radial  distribution  function  for  3D  moderately  dense- 
packed  configurations  of  90  spheres(discrete  points),  compared  with  the  P-Y  radial 
distribution  function(solid  curve) 


Chapter  4 


Experimental  Investigation 


The  following  sections  describe  the  triaxial  compression  experiment  employ¬ 
ing  steel  baUs  as  conductive  granular  particles. 


4.1  Equipment 

4.1.1  Compression  Tester 

The  current  experiments  employ  a  commercial  compression  tester  as  loading 
frame  for  the  test  cell.  The  911  MTT-02/10  compression  tester,  manufactured  by 
Comten  Industries  Inc.,  is  a  motorized  device  with  a  digital  interface  and  a  variable 
speed  drive.  It  consists  of  two  parts,  the  m2un  unit  and  a  monitor/controller.  The 
loading  force  is  measured  and  displayed  on  the  monitor /controller.  The  displacement 
of  a  specimen  is  measured  with  a  separate  500  DC-E  LVDT  linear  displacement 
transducer  manufactured  by  Lucas-Schaevitz. 

4.1.2  Triaxial  Cell 

The  triaxial  cell,  a  redesigned  version  of  standard  commercial  cell  with  ad¬ 
ditional  provision  for  conductivity  measurement,  is  schematically  illustrated  in  Fig¬ 
ure  4.1.  Furthermore,  the  cell  dimensions  have  been  reduced  in  order  to  fit  into  the 
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compression  tester,  and  a  bellows  is  added  to  the  top  cap  to  prevent  escape  of  con¬ 
fining  fluid  and  to  eliminate  the  friction  between  piston  (shaft)  and  its  bushing,  a 
problem  that  is  inevitable  in  the  conventional  design  of  the  triaxial  cells.  The  cell 
is  made  of  stainless  steel,  with  a  transparent  plexiglass  cylindrical  body.  The  base 
adaptor,  with  a  porous  steel  screen  installed  on  the  top,  has  a  channel  for  connection 
through  the  cell  base  to  a  vacuum  system,  for  initial  specimen  preparation,  or  to 
the  atmosphere  for  drained  tests.  The  loading  cap  h2is  also  a  porous  steel  screen  at 
the  bottom  and  an  eccentric,  small-diameter  hole  on  the  top  to  rele£ise  air  during 
the  preparation  of  saturated  specimens  (not  used  in  the  current  experiments).  This 
small- diameter  hole  is  always  blocked  by  the  threaded  plug  during  experiments  to 
prevent  compressed  air  from  escaping  through  the  specimen.  In  the  center  of  the 
loading  cap,  there  is  a  conical  recess,  into  which  the  end  of  loading  shaft  seats.  In  the 
top  cap,  there  is  an  outlet  from  which  the  compressed  air  enters  the  cell  chamber  to 
form  a  confining  pressure  around  the  specimen.  The  LVDT  device  is  mounted  on  the 
top  cap  to  measure  the  end  displacement  of  the  specimen.  The  external  connection 
to  measure  the  resistance  through  the  granular  specimen  is  also  shown  in  Fig.  4.1. 

4.1.3  Digital  Image  Processing  System 

A  PC-based  digital  image  acquisition/processing  system,  used  in  a  previous 
study  [140],  is  employed  for  data  acquisition.  The  hardware  consists  of  a  Hitachi 
KVC-150  video  camera,  a  Toshiba  M-6100  VCR,  a  Sony  PVM-1271Q  monitor,  a 
IBM-PC  compatible  and  an  embedded  PCVISIONplus  PFGPLUS-512-3-U-XT/AT 
frame  grabber,  manufactured  by  Image  Technology  Inc..  The  software  includes  a 
Werner-Frei  Associates  Image  Lab  and  Imagetool  programs. 

The  basic  components  of  the  system  and  their  mutual  relation  are  depicted 
in  Figure  4.2. 
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PCVISIONplus  Frame  Grabber 


Figure  4.2:  Digital  acquisition/processing  system,  (after  Zhuang  1991). 

4.2  Materials  and  Specimens 

The  test  material  used  as  particles  in  the  triaxial  compression  experiments 
is  440-C  stainless  steel  balls  provided  by  Thomson  Precision  Ball  Company,  Inc.,  the 
physical  properties  being  listed  in  Table  4.1  in  which  the  first  four  properties  are 
provided  by  the  manufacturer,  the  fifth  is  based  on  our  own  measurement,  and  the 
last  one  is  found  from  CRC  Handbook  of  the  tables  for  Applied  Engineering  Science 
[21],  An  interparticle  friction  coefficient  /i  =  0.15  was  meaisured  by  gluing  the  steel 
balls  to  a  plate  and  observing  the  critical  angle  of  shding  down  a  second  inchned 
stainless  steel  plate.  The  granular  assemblages  consisted  of  a  randomly  packed  beds 
of  approximate  3400  of  the  above  steel  balls  having  a  3.81cm  diameter  and  a  4.9cm 
height.  The  initial  packing  densities  were  around  0.6,  which  is  close  to  that  of  a 
random  dense  packing  of  spheres.  The  ambient  confining  pressure  is  kept  constant  at 
QA^kglcm?  throughout  the  course  of  the  deformation. 
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Diameter  (cm) 

0.278  ±  0.000064 

Density  {gfcm^) 

7.667 

Elastic  modulus  (kg/crn^) 

2.039  X  10® 

Poisson  ratio 

0.29 

Friction  coeflBcient 

0.15 

Electrical  conductivity  (ohm  •  cm)“^ 

16670.0 

Table  4.1:  The  physical  properties  of  the  steel  balls  employed  in  the  experiments. 

4.3  Specimen  Preparation  and  Experimental  Pro¬ 
cedure 

First,  the  latex  membrane  is  attached  to  the  base  adaptor,  then  stretched 
to  the  top  of  the  split-cylinder,  two-piece  mold  to  form  a  cylindrical  space.  The  steel 
balls  are  poured  into  the  mold  space  in  five  portions,  each  is  followed  by  rod-stirring, 
a  procedure  known  as  “rodding”  in  the  soil-mechanics  community,  in  order  to  achieve 
a  consistent  initial  packing  condition.  The  loading  cap  is  then  placed  on  the  specimen 
and  held  by  the  membrane.  A  vacuum  is  then  applied  to  withdraw  the  air  from  the 
specimen,  which  makes  the  specimen  rigid  under  the  atmospheric  pressure.  After  the 
two-piece  mold  is  removed,  the  plexiglass  cylindrical  body  together  with  the  top  cap 
is  installed  on  the  cell  base  and  tightened  with  three  tie  bars.  The  cell  is  then  placed 
between  two  platens  of  the  compression  tester.  The  LVDT  transducer  is  mounted 
on  the  top  cap  of  the  cell.  Laboratory  compressed  air  is  used  to  fill  an  approximate 
O.OSm^  tank  serving  as  air  reservoir,  which  maintains  a  constant  confining  pressure 
during  the  experiments.  Simultaneously,  a  confining  pressure  is  created  inside  the 
triaxial  cell  that  is  always  connected  to  the  jur  reservior.  The  compressed  air  is  shut 
off  when  the  desired  pressure  p  is  reached.  The  vacuum  is  then  turned  off,  and  a 
valve  is  opened  to  vent  the  air  in  the  specimen  to  accommodate  drained  tests.  The 
specimen  preparation  is  now  complete. 

Following  the  above  specimen  preparation,  the  compression  tester  is  turned 
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on.  The  wires  are  connected  between  the  cell  and  one  electrical  multimeter,  which 
is  utilized  to  measure  the  electrical  resistance  through  the  specimen.  The  LVDT 
is  connected  to  a  second  multimeter,  which  measures  the  voltage  output  from  tlie 
LVDT,  thereafter  being  converted  to  the  axial  displacement.  The  video  camera  and 
VCR  are  set  up  to  record  the  display  readings  of  the  monitor /controller  and  two 
multimeters.  Everything  is  then  ready  for  the  experiment. 

VCR  recording  and  the  compression  tester  are  then  started.  The  loading 
speed  is  slowly  increased  from  zero  to  approximately  O.lcm/min.  This  loading  process 
is  terminated  at  20%  axial  compressional  strain. 


4.4  Data  Analysis 


With  the  help  of  the  image  processing  system,  about  30  select  frames  are 
extracted  from  the  recorded  video  image  of  the  entire  experiment  at  different  stages, 
and  the  force,  resistance  and  voltage  readings  are  read  off  for  later  data  analysis. 

The  measured  force,  corrected  for  the  bellows  spring  force,  gives  us  the  axial 
force  exerted  on  the  specimen  by  the  shaft.  Then  we  compute  the  vertical  stress  over 
the  specimen  according  to: 


O’!  =  p  + 


(4.1) 


where  p  is  the  confining  pressure,  F  is  the  axial  force  exerted  on  the  specimen  by  the 
shaft,  and  Acap  is  the  area  of  the  loading  cap.  The  horizontal  stresses  <72  and  <73  are 
equal  to  the  confining  pressure  p. 

According  to  the  definition,  we  obtain 


a:*  = 


L 

RtAg 


(4.2) 


where  Rt  is  the  total  resistance  through  the  specimen,  K*  is  the  effective  conductivity 
of  the  medium,  L  is  the  specimen  length  and  A,  is  the  cross-section  area  of  the 
specimen. 
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4.5  Measurements  of  Contact  Resistance 

Faced  with  the  extremely  large  discrepancy  between  the  experimental  mea¬ 
surements  and  the  results  of  the  initial  numerical  simulations  utilizing  eq.  2.23  to 
compute  the  contact  resistance,  we  were  lead  to  identify  the  underlying  cause.  As 
a  matter  of  fact,  there  has  been  ongoing  research  into  electric  contact  resistance  for 
the  past  several  decades  [64,60,61,6,125,68,83,24,119,25],  since  almost  all  electrical 
or  electronic  equipment  contains  numerous  contacts  through  which  the  current  and 
voltage  signals  are  transmitted,  and  the  failure  of  even  a  single  contact  can  result 
in  a  complete  system  failure.  From  the  above  literature  survey,  one  finds  that  the 
resistance  at  a  real  contact  may  be  much  higher  owing  to  the  various  reaisons  outlined 
in  the  section  2.3.2.  This  ultimately  results  in  lower  effective  conductivity  of  granular 
medium  than  that  predicted  by  the  theory. 

Based  on  this  understanding,  we  conducted  an  experiment  to  measure  the 
electrical  contact  resistance  in  a  single  column  of  steel  balls  as  a  function  of  axial 
load.  The  e.xperimental  set-up  is  schematically  depicted  in  figure  4.3.  A  stack  of  balls 
is  confined  within  a  reinforced  quartz  glass  capillary  tube,  with  an  inner  diameter 
0.29cm,  slightly  larger  than  that  of  the  steel  balls,  and  loaded  with  dead  weight 
from  the  top  via  a  piston.  The  resistance  is  meeisured  with  a  multimeter  at  various 
loads.  To  avoid  all  contacts  between  the  balls  and  flat  surfaces,  the  top  and  bottom 
steel  balls  in  the  column,  consisting  of  five  balls,  are  soldered  to  the  piston  and  the 
flat  baise,  respectively.  The  relationship  between  resistance  and  load,  the  average 
result  of  seven  experiments,  is  plotted  in  figure  4.4  and  compared  with  that  predicted 
by  Eq.  2.23.  The  overall  scatter  is  within  20%  of  the  average  values.  The  large 
scatter  is  probably  ascribed  to  the  nonuniformity  of  the  ball  surfaces  and  associated 
contamination  films.  One  can  see  a  very  strong  dependence  of  resistance  on  the  load, 
with  a  slope  of  approximately  2.4  compared  with  1/3  given  by  Hertzian  theory.  This 
observation  suggests  that  an  oxide  film  on  the  ball  surface  ruptures  and  deforms 
plastically  under  the  applied  load. 
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Figure  4.3:  Experimental  set-up  used  to  measure  the  electrical  contact  resistance. 

The  load-resistance  curve  depicted  in  Fig.  4.4  represents  the  case  of  loading. 
Experimental  observations  revealed  that  the  load-resistance  curve  for  unloading  is 
lower  than  that  of  loading.  This  hysteresis  is  probably  due  to  the  plastic  deformation 
of  asperities  and  irreversible  rupture  of  the  superficial  oxide  film  during  loading. 
While  the  contact  between  two  balls  can  experience  loading  and  unloading  during  the 
deformation  of  granular  media,  the  localized  contact  regions  on  the  ball  surfaces  keep 
changing  due  to  the  relative  movement  (rolling  and  sliding)  between  balls.  Therefore, 
it  is  plausible  to  employ  the  normal  load-contact  resistance  relation  of  Fig.  4.4  for 
loading  for  purposes  of  numerical  simulation. 


of  the 


Chapter  5 


Simulation  Results 


A  series  of  numerical  simulations  on  idealized  granular  assemblages  have 
been  conducted  to  investigate  the  effects  of  microscopic  and  microstructural  properties 
of  the  constituent  particles  and  their  assemblages,  including  (Coulomb)  interpar  ‘  Icle 
friction,  nonlinear  contact  mechanics  and  initial  packing  density,  on  the  mechani¬ 
cal  behavior.  Of  particular  interest  is  the  Reynolds  dilatancy,  shear  strength  and 
the  evolution  of  granular  microstructure  of  idealized  granular  aissemblages  subject  to 
constant  mean  confining  pressure.  Numerical  simulations  of  the  triaxial  compression 
test  have  been  conducted  to  simulate  the  effects  of  the  initial  density  on  the  mechan¬ 
ical  behavior  as  well  as  the  scalar  transport  properties.  The  mechanical  behavior  is 
discussed  first  in  this  chapter,  while  transport  properties  will  be  covered  separately 
in  the  following  chapter. 


5.1  Interparticle  Friction 

This  study  involves  both  2D  and  3D  mono-  and  poly-disperse  granular  as¬ 
semblages  subject  to  simple  shear  deformation  under  constant  mean  confining  pres¬ 
sure.  The  2D  assemblages  consist  of  132  disks  initially  packed  to  random  dense 
packing  with  about  0.82  area  fraction.  On  the  other  hand,  the  3D  granular  assem¬ 
blages  contain  48  spheres  initially  packed  to  an  approximate  dense  random  packing 
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with  0.60  volume  fraction.  The  nondimensional  radii  of  the  particles  are  equal  to 
1.0  for  monodisperse  systems,  and  equal  to  0.8,  1.0,  and  1.25,  respectively,  in  the 
polydisperse  systems,  with  approximately  same  total  volume  of  particles  of  the  three 
different  sizes.  Thus,  we  use  64  disks  of  radius  0.8,  41  disks  of  radius  1.0  and  27 
disks  of  radius  1.25  for  2D  systems,  and  choose  27  spheres  of  radius  0.8,  14  spheres 
of  radius  1.0  and  7  spheres  of  radius  1.25  for  3D  systems.  Normal  and  shear  contact 
stiffnesses  are  taken  to  be  kn  =  10  and  kt  =  0.8,  respectively. 

By  a  scaling  based  on  contact  stiffness  and  particle  radius,  one  can  specify 
an  externally  imposed  nondimensional  pressure  (=p,i2r/fc„,,  where  pr,  kn,  ®nd  Rr 
denote  the  real  confining  pressure,  normal  stiffness  and  particle  radius.),  under  which 
interparticle  overlap(proportional  to  normal  force)  will  not  exceed  0.1%  of  particle 
radius  throughout  the  simulations  since  we  are  primarily  concerned  with  nearly  rigid 
particles  [14].  This  pressure  is  found  to  be  6.0  x  10“®  for  2D  experiments  and  4.0  x 
10”®  for  3D  experiments,  and  is  maintained  during  the  initial  packing  process  and 
subsequent  shearing.  Both  the  2D  or  3D  test  assemblages  are  subjected  to  simple 
shear  up  to  20%  total  strain  with  different  interparticle  friction  coefficients  under 
otherwise  identical  conditions. 

To  further  clarify  the  influence  of  particle  friction  on  Reynolds  dilatancy  of 
randomly  dense-packed  granular  a  ^  jmblages,  we  carried  out  several  simulations  on 
both  2D  and  3D,  and  mono-  as  well  as  poly-disperse  idealized  granular  assemblages 
with  /i=0.0,  0.3,  and  0.5,  respectively.  The  following  conclusions  can  be  drawn  from 
the  results,  presented  in  Figure.  5.1,  5.2,  5.3,and  5.4:  the  dilatancy  increase  with 
increasing  magnitude  of  /x,  which  is  in  agreement  with  previous  results  [39,18,33], 
including  the  results  for  polydisperse  (random-packed)  cases  found  by  Bashir  and 
Goddard[14].  However,  this  finding  is  contrary  to  Reynolds’  original  hypothesis  on 
the  random  dense  packing  of  granular  assemblages,  as  interpreted.  The  stress  ratio 
(o-i  —  <T3)fp,  where  ci,  (Ts  and  p  are  major,  minor  and  mean  stress,  also  increases  with 
increasing  magnitude  of  p.  Polydispersity  is  found  to  have  a  noticeable  effect  on  the 
mechanical  behavior  of  granular  assemblages. 


50 


From  Fig.  5.5  and  5.6,  one  finds  that  the  coordination  number  decreased 
drastically  in  the  start  of  shearing,  usually  within  1%  of  shear  strain,  which  indicates 
that  a  significant  particle  rearrangement  took  place  early  in  the  shearing  deformation 
[79,95,96].  The  coordination  number  then  fluctuated  about  a  roughly  constant  level 
throughout  the  subsequent  deformation  [18].  A  higher  interparticle  friction  generally 
results  in  a  lower  coordination  number  after  granular  assemblages  yield.  Although 
the  computed  coordination  number  varies  with  shear  strain  and  interparticle  friction, 
it  is  fotmd  to  be,  in  both  2D  and  3D  cases,  always  greater  than  the  geometric  critical 
coordination  number  Zc  «  (=  2  for  d  =  2  and  1.5  for  d  =  3,  here,  d,  >  1,  denotes 

the  number  of  dimensions  [116]).  at  which  there  exists  an  infinite  cluster  of  bonds 
across  a  medium  according  to  percolation  theory  [116,75],  is  a  dimensional  invariant 
insensitive  to  the  details  of  the  lattice  studied.  The  geometric  percolation  threshold  pc 
is  shown  to  be  about  0.347  for  2D  triangular  lattice  and  0.119  for  FCC  lattice  [116], 
while  the  elastic  (central-force,  omitted  hereafter)  bond  percolation  threshold  pcen> 
which  would  associate  with  solid  behavior  of  bond  networks  at  small  strain,  is  found 
to  roughly  equal  to  0.58  for  2D  triangular  lattice  and  0.42  for  FCC  lattice  [49].  Here, 
we  define  a  ratio  ZjZmax  representing  the  fraction  of  active  “bonds”  in  the  network 
of  particle  contacts  compared  with  the  coordination  number  of  the  densest  possible 
systems,  where  Z  being  the  coordination  number  of  a  system  and  Z^ax  maximum 
possible  coordination  number,  for  instance,  6  for  2D  triangular  lattice  and  12  for  3D 
FCC  lattice.  In  2D  monodisperse  case,  from  Fig.  5.5,  one  finds  a  ratio  0.6  at  initial 
stage,  slightly  larger  than  the  elastic  bond  percolation  threshold,  and  a  range  from 
0.4  to  0.5,  dependent  on  interparticle  friction,  after  early  shearing  deformation,  which 
is  between  the  geometric  and  elastic  percolation  threshold.  Similar  results  are  found 
in  3D  monodisperse  case.  In  a  previous  2D  work,  however,  Bashir  and  Goddard  [14] 
found  the  ratio  very  close  to  the  geometric  percolation  threshold  after  initial  3%  shear 
strain,  which  is  much  smaller  than  the  elastic  percolation  threshold. 

Based  on  detailed  microscopic  observations,  we  find  that:  (1)  granular  mi¬ 
crostructure  evolves  such  that  contact  normals  concentrated  in  the  direction  of  major 
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Figure  5.1:  Effects  of  interpaxticle  friction  coefficients  on  dilatancy  of  2D  assemblages 
subjected  to  simple  shearing  deformation. 

principal  stress  during  the  shearing  deformation  [124];  (2)  the  granular  assemblage  is 
composed  of  two  types  of  region:  a  major  skeleton  composed  of  heavily  stressed  chains 
of  particles  and  less  stressed  regions  surrounding  this  skeleton,  with  most  of  contact 
breaking  and  making  occurring  within  the  latter  region  and  with  the  skeleton  remain¬ 
ing  relatively  unaltered  for  a  small  incremental  deformations;  (3)  particle  rolling  is 
major  deformation  mechanism,  especiadly  when  interparticle  friction  is  large[102]. 

5.2  Nonlinear  Contact  Mechanics 

As  we  discussed  in  section  2.2,  the  contact  stiffness  is  generally  a  function 
of  load,  often  well  represented  by  the  power  law 

K  =  Cfi  (5.1) 


where  C  is  a  material  constant,  fn  normal  load,  and  for  example,  the  exponent  takes 
on  values  A  =  1/3  for  Hertzian  elastic  spheres  and  A  =  1/2  for  a  conical  tip  pressed 
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Figtire  5.2:  Effects  of  interparticle  friction  coefficients  on  dilatancy  of  3D  assemblies 
subjected  to  simple  shearing  deformation. 


Figure  5.3:  Effects  of  interparticle  friction  coefficients  on  shear  strength  of  2D  assem¬ 
blages  subjected  to  simple  shearing  deformation. 


Figure  5.4:  Effects  of  interparticle  friction  coefficients  on  shear  strength  of  3D  assem¬ 
blages  subjected  to  simple  shearing  deformation. 


Figure  5.5:  Effects  of  interparticle  friction  coefficients  on  average  coordination  number 
of  2D  assemblages  subjected  to  simple  shearing  deformation. 
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Figure  5.6:  Effects  of  interparticle  friction  coefficients  on  average  coordination  number 
of  3D  assemblages  subjected  to  simple  shearing  deformation. 

against  a  plane  [57].  Accordingly,  the  relationship  3.3  between  increments  of  relative 
displacement  and  contact  force  is  nonlinear. 

Since  the  linear  model  A  =  0  not  only  offers  simplicity  but  may  be  able  to 
provide  qualitatively  valid  insights  into  the  link  between  micromechanical  properties 
and  macroscopic  behavior  [17,18,39,40,41,128]  and,  since  we  are  mainly  interested 
in  ideal  rigid  limit  [14],  most  of  our  simulations  have  been  carried  out  for  linear 
contacts.  However,  we  felt  that  it  is  important  to  check  its  validity  and  the  effects  of 
nonlinearity. 

In  the  present  simulations  we  employ  a  monodisperse  system  with  48  spheres 
packed  to  initial  density  <f>  =  0.60  and  interparticle  friction  (i  =  0.15.  The  system 
is  subjected  to  the  triaxial  cor"  'ssion  under  a  constant  mean  confining  pressure 
p  =  4  X  10“*.  A  values  are  selected  as  0.0,  1/3,  1/2,  and  1,  with  0.0  representing  the 
linear  contact  and  1.0  representing  extreme  nonlinearity.  The  tangential  stiffness  kt 
is  simply  taken  to  be  O.Skn. 

From  Fig.  5.7,  one  sees  that  contact  nonlinearity  has  no  apparent  influence 
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Figure  5.7:  Effects  of  non-linear  contact  on  dilatancy  of  3D  assemblages  subjected  to 
I  triaxial  compression. 
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Figure  5.8:  Effects  of  non-linear  contact  on  average  coordination  number  of  3D  as- 
I  semblages  subjected  to  triaxial  compression. 
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Figure  5.9:  Effects  of  non-linear  contact  on  average  contact  normal  force  /„  in  3D 
assemblages  subjected  to  triaxial  compression. 


Figure  5.10:  Effects  of  non-linear  contact  on  shear  strength  of  3D  assemblages  sub¬ 
jected  to  triaxial  compression. 
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on  Reynolds  dilatancy  in  the  small-  to  intermediate-strain  region,  but  some  effects 
are  observed  at  higher  axial  strains.  Moreover,  Fig.  5.8  indicates  that  the  average 
number  of  particle  contacts  is,  roughly  speaking,  not  affected.  It  is  expected  that  the 
nonlinear  contact  law  will  tend  to  make  the  strong  contacts  (in  terms  of  contact  force) 
stronger,  and  the  weak  contacts  weaker,  therefore,  changing  the  force  distribution. 
Although  we  do  indeed  observe  variations  in  the  force  distributions,  the  different 
degrees  of  nonlinearity  yield  only  a  small  deviation  in  the  average  normal  force,  for 
the  cases  A  =  0.0,  1/3  and  1/2  (see  Fig.  5.9).  Finally,  from  the  results  in  Fig.  5.10, 
one  sees  that  the  strength  of  the  idealized  granular  assemblages  tends  to  increase 
somewhat  with  the  increase  of  nonlinearity.  The  small  effects  of  contact  nonlinearity 
on  the  mechanical  behavior  indicate  that  a  complicated  nonlinear  contact  law  is  less 
significant  in  modeling  the  mechanical  behavior  of  granular  materials. 


5.3  Effects  of  Initial  Specimen  Density 

To  simulate  the  effects  of  the  initial  void  ratio  or  density  on  the  behavior 
of  mechanical  as  well  as  the  scalar  transport  properties,  we  have  generated  three 
random  monodisperse  packings  of  100  spheres  with  different  initial  densities,  0.52, 
0.56  and  0.60,  respectively.  The  interparticle  friction  coefficient  /x  is  taken  to  be  0.15. 
All  three  packings  are  subjected  to  the  triaxial  compression  deformation  under  the 
same  nondimensional  confining  pressure  po  =  4  x  10“®  in  the  directions  normal  to 
compressional  axis. 

Fig.  5.11  and  5.12  indicate  that,  for  initieilly  loose  systems  such  as  those 
with  <l>  =  0.52  and  0.56,  densification  or  negative  dilatancy  occurs  initially  and  per¬ 
sists  throughout  the  deformation.  The  potential  for  densification  increases  with  the 
decrease  of  the  initial  density  [91].  On  the  other  hand,  the  initially  dense  system, 
with  <f>  =  0.60,  experiences  positive  dilatancy  from  the  very  beginning  of  the  defor¬ 
mation.  Nevertheless,  the  densities  of  three  systems,  either  contracting  or  expanding, 
tend  to  approach  the  same  critical  value  asymptotically.  Fig.  5.13  shows  that,  for 
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loose  systems,  the  shear  strength  increases  monotonic2Jly.  However,  for  the  initially 
dense  system,  its  strength  increases  first  until  a  peak  value  is  reached,  then  decreases. 
Similar  observations  in  real  triaxial  compression  tests  are  to  be  discussed  in  follow¬ 
ing  chapter.  Again,  both  the  loose  and  dense  systems  seem  to  possess  an  identical 
ultimate  strength  after  being  subjects!  to  a  large  deformation  [91].  The  initial  coor¬ 
dination  number  increases  as  the  density  rises.  Upon  deformation,  the  dense  system 
initially  experiences  a  significant  loss  in  the  number  of  contacts,  whereas  the  loose 
system  gains  contacts.  However,  both  systems  approach  approximately  the  same 
coordination  number  at  roughly  1%  axial  compressional  strain,  and  then  maintain 
fluctuating  but  slightly  increasing  values  (see  Fig.  5.14).  Next,  I  shall  attempt  to 
elucidate  the  observed  behavior. 

In  a  loose  system,  the  initial  number  of  interparticle  contacts  and  coordi¬ 
nation  number  are  low,  and  just  exceeds  slightly  the  elastic-percolation  threshold,  at 
which  there  just  begin  to  exist  sample-spanning  chains  of  particles  capable  of  support¬ 
ing  an  ambient  confining  stress  (Goddard  1990).  However,  due  t-'  lack  of  sufficient 
contact  force  from  neighboring  particles,  such  chains  are  highly  unstable  to  (Euler) 
buckling.  Under  these  circumstances,  one  can  anticipate  that  a  given  pair  of  adjacent 
particles  in  such  a  chain  will  accommodate  axial  compression  with  a  small  rotation 
normal  to  their  line  of  centers  until  such  rotation  is  hindered  by  lateral  contact  with 
neighboring  particles.  By  means  of  this  process,  a  given  particle  chain  will  generally 
undergo  a  kind  of  lateral  ‘branching’  until  it  becomes  capable  of  supporting  increased 
axial  compressive  stress  (Goddard  1990).  Therefore,  the  overall  granular  structure  is 
less  stable  and  more  likely  to  collapse  to  a  more  stable,  denser  system  upon  deforma¬ 
tion  and  to  generate  load-bearing  capability.  Such  capability  is  further  enhanced  as 
the  system  gets  denser.  On  the  other  hand,  the  dense  system  wiU  have  to  expand  in 
order  to  deform,  hence  loses  contacts  initially.  Owing  to  the  volume  exoansion  against 
the  ambient  confinement,  the  system  exhibits  shear  strength,  but  further  dilatancy 
reduces  the  system  density  and  decrezises  the  stability  of  the  granular  chain  structure 
and  its  ability  to  support  the  external  loads.  This  explains  the  after-peak  strength 
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Figure  5.12:  Effects  of  initial  density  on  density  evolutions  of  3D  assemblages  sub¬ 
jected  to  triaxial  compression. 


Figure  5.13:  Effects  of  initial  density  on  shear  strength  of  3D  assemblages  subjected 
to  triaxial  compression. 


Chapter  6 


Electrical  Conductivity 


6.1  Numerical  Simulations 

For  the  purpose  of  data  validation,  the  computer  code  is  modified  specifically 
to  simulate  a  triaxial  compression  experiment  over  an  idealized  granular  assemblage 
with  the  physical  properties  of  steel  balls,  listed  in  table  4.1 ,  and  experimental  loading 
conditions  as  input  parameters.  The  mechanical  behavior  as  well  as  the  electrical 
transport  properties  of  the  systems  are  investigated  simultaneously. 

The  idealized  system,  as  schematically  depicted  in  figure  6.1,  contains  100 
equal  size  spherical  particles  in  a  periodic  cubic  cell.  The  system,  being  confined  in 
X  and  Y  directions  with  a  constant  pressure  po,  is  compressed  in  Z  direction.  To 
simulate  the  effects  of  the  initial  density  on  the  mechanical  eis  well  as  scalar  transport 
properties,  we  have  generated  three  isotropic  random  packings  of  sphere  2issemblages 
with  initial  densities  being  0.52,  0.56  and  0.60,  respectively.  The  interparticle  friction 
coefficient  p  is  taken  to  be  0.15,  corresponding  to  the  measurements  on  the  real 
‘dirty’  steel  balls  (whose  definition  is  to  be  provided  in  Section  6.2).  The  computed 
dimensionless  normal  force  is  converted  to  real  force  by  the  scaling  factor  based  on  the 
particle  radius  and  ambient  confining  pressure  so  that  the  contact  resistance  between 
two  particles  can  be  calculated  with  the  experimentally  determined  load-resistance 
relation  given  in  figure  4.4,  also  corresponding  to  the  ‘dirty’  steel  balls.  The  tangential 
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Figure  6.2:  Relation  between  mean  effective  conductivity  and  specimen  density. 

force  is  not  considered  in  evaluating  the  contact  resistance,  since,  according  to  our 
experimental  observations,  the  tangential  force  has  no  apparent  effert  on  the  contact 
resistance.  The  results  reported  in  this  section  correspond  to  the  physical  properties 
of  ‘dirty’  steel  balls. 

It  is  generally  true  that  a  denser  system  should  possess  higher  conductivity  or 
lower  resistivity  owing  to  the  fact  that  medium  is  better  connected  so  Jis  to  offer  more 
paths  or  branches  for  current  to  pjiss  through.  Fig.  6.2,  where  the  mean  effective  (the 
word  ‘effective’  and  corresponding  superscript  ‘*’  is  henceforth  omitted)  conductivity 
is  the  average  of  conductivities  Kn,  Kyy  and  Kzz  in  X,  Y  and  Z  directions,  reveals 
this  trend  although  more  fluctuation  are  observed  in  high  density  region. 

Next,  we  shall  consider  the  changes  of  the  conductivity  Kzz  in  the  compres- 
sional  (or  Z)  direction  during  the  course  of  deformation.  From  Fig.  6.3  and  5.13,  it  is 
not  difficult  to  see  the  similarity  between  the  behavior  of  shear  strength  and  conduc¬ 
tivity,  although  there  are  more  fluctuations  in  the  conductivity  than  in  strength.  For 
the  initially  dense  system,  with  <{>  =  0.60,  the  conductivity  increases  first  to  a  peak 
within  the  first  2%  of  axial  strain,  then  fluctuates  wildly  with  a  decreasing  trend. 
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Figure  6.3:  Dependence  of  conductivity  on  axial  strain  and  initial  density. 

This  behavior  can  be  correlated  to  the  mechanical  behavior,  particularly  deformation 
mechanism  of  the  granular  assemblages. 

Upon  the  compression,  load-bearing  chains  are  built  up  gradually  in  the 
Z  direction  in  order  to  sustain  stress  anisotropy,  which  create  easier  paths  for  cur¬ 
rent  to  pass  through.  Therefore,  in  early  stage  of  deformation,  the  conductivity 
increases  steadily  despite  the  fact  that  system  experiences  a  loss  in  total  number  of 
contacts  (see  Fig.  5.14),  note  that  the  loss  of  contacts  happens  mainly  on  X-Y  plane 
or  the  direction  of  minor  principal  stress.  Owing  to  the  dilatancy  of  the  system  (see 
Fig.  5.11),  these  granular  chains  become  less  stable  progressively.  When  the  system 
is  further  expanded,  these  major  load-bearing  chains  finally  buckle.  The  branching- 
out  of  chains  diverts  current  from  the  preferred  direction.  Therefore,  one  observes 
the  after-peak  decrease  in  conductivity.  The  subsequent  built-up  and  buckling  of 
new  heavily-stressed  chains  is  primary  cause  of  the  fluctuations  in  conductivity.  On 
the  other  hand,  the  densification  in  the  loose  systems  during  the  deformation  tends 
to  stablize  these  progressively  loaded  chains,  thus  resulting  in  a  steady  increase  of 
conductivity. 


e. 


Figure  6.4;  Evolution  of  principal  value  ratios  for  stress,  fabric  and  conductivity 
=  0.60). 

The  principal  value  ratio,  the  ratio  of  major  to  minor  principal  value  of 
tensors  such  as  stress  or  fabric,  is  often  used  to  characterize  the  state  of  anisotropy. 
Following  is  an  attempt  to  correlate  the  evolution  of  the  anisotropy  of  three  tensors, 
namely  stress  <r,  fabric  N  and  conductivity  K  tensors.  At  the  start  of  compression, 
the  anisotropy  of  the  assemblages  is  induced  progressively.  After  reaching  a  peak, 
it  remains  relatively  unchanged  throughout  the  deformation.  From  Fig.  6.4,  one  see 
fabric  is  not  very  sensitive  to  the  change  in  stress  state,  while  conductivity  is  more 
sensitive. 

Fig.  6.5  and  6.6  clearly  suggest  strong  correlations  between  conductivity, 
stress  and  fabric  tensors.  Note  that  the  lines  of  linear  regression  almost  pass  through 
the  point  (1,1)  in  the  plots,  which  represents  the  isotropic  relation  between  these 
tensors. 

The  results  shown  in  Fig.  6.7  indicate  that  the  mean  field  theory  of  Batch¬ 
elor  and  O’Brien  underpredicts  the  effective  conductivity  slightly  but  can  be  used  to 
understand  the  qualitative  scalar  transport  behavior  of  granular  materials. 


Figure  6.5:  Relations  between  principal  value  ratios  of  conductivity  and  stress 
(slope=2.39). 


Figure  6.6:  Relations  between  principal  value  ratios  of  conductivity  and  fabric 
(slope=5.23). 
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Figure  6.7:  Comparison  between  the  exact  solution  and  mean  field  theory  in  predict¬ 
ing  the  effective  conductivity  of  idealized  granular  assemblages  {<i>  =  0.60). 

6.2  Experimental  Verifications 

Three  triaxial  compression  tests  have  been  carried  out,  one  with  ‘clean’  balls 
and  two  with  ‘dirty’  balls,  at  approximate  <f>  =  0.60  initial  densities.  Here,  clean  balls 
refer  to  those,  whose  protective  oil  film  wm  just  removed  and  cleaned  by  Aceton  sol¬ 
vent,  and  are  expected  to  have  no  much  contamination  film  on  the  ball  surface  (or  low 
contact  resistance,  however,  the  normal  load-contact  resistance  relation  was  not  ex¬ 
perimentally  determined  immediately  after  cleaning  process  unfortunately),  whereas 
the  dirty  balls  refer  to  those,  which  have  been  exposed  to  normal  laboratory  environ¬ 
ment  for  approximately  four  months  after  being  cleaned  with  Aceton,  thereby,  possess 
thicker  insulating  contamination  film  (or  high  contact  resistance),  and  for  which  the 
load-contact  resistance  relation  (Fig.  4.4)  and  interparticle  friction  coefficient  have 
been  measured.  The  experimental  results  with  clean  balls  are  given  here  solely  for 
the  comparison  against  the  experimental  results  with  dirty  balls.  It  is  not  intended 
to  compare  the  experimental  results  with  clean  balls  to  the  results  of  the  numerical 
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Figure  6.8:  Comparison  of  shear  strength  between  numerical  simulation  and  experi¬ 
ments  with  dirty  balls  (<f>  =  0.60). 

simulations,  since  no  numerical  results  for  clean  balls  are  available. 

Plots  Fig.  6.8  and  6.9  indicate  that  the  results  of  numerical  simulation  and 
physical  experiments  with  dirty  steel  balls  are  in  qualitative  agreement,  and  that  the 
simulation  is  capable  of  predicting  the  mechanical  and  scalar  transport  properties 
of  granular  assemblages.  Comparison  of  the  experimental  results  between  clean  and 
dirty  balls  in  Fig.  6.9  also  reveals  that  the  individual  contact  resistance  can  drastically 
affect  the  conductivity  of  the  medium. 
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Figure  6.9:  Comparison  of  electrical  conductivity  between  numerical  simulation  and 
experiments  with  dirty  balls  =  0.60). 
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Chapter  7 

Conclusions  and 


Recommendations 


We  have  developed  a  new  version  of  a  quasi-static  simulation  of  the  me¬ 
chanics  and  conductivity  of  sphere  assemblages  by  introducing  several  new  numerical 
techniques,  including  a  relaxation  method  which  is  shown  to  be  a  useful  tool  to  over¬ 
come  a  singularity  in  the  quasi- linear  system  of  equilibrium  equations.  The  computer 
code  is  versatile  enough  to  allows  one  to  simulate  any  deformation  history  and  to 
study  both  mechanical  and  scalar  transport  properties  of  idealized  granular  assem¬ 
blages  simultaneously. 

The  results  of  the  present  investigation  show  that:  (1)  interpartide  friction 
has  great  influence  on  Reynolds  dilatancy  for  random  dense-packed  mono-  as  well  as 
poly-disperse  granular  assemblages,  a  result  contrary  to  Reynolds’  original  hypothesis; 

(2)  the  use  of  linear  contact  mechanics  is  justified  near  the  ideal  rigid-particle  limit; 

(3)  scalar  transport  properties  such  &s  electrical  conductivity  can  be  employed  as  a 
good  indicator  of  the  stress  anisotropy  and  microstructural  particle-contact  topology, 

(4)  the  comparison  between  numerical  results  and  experimental  findings  reveals  that 
the  numerical  model  is  able  to  qualitatively  predict  the  mechanical  as  well  as  scalar 
transport  properties  of  idealized  granular  assemblages;  (5)  the  contact  resistance  be¬ 
tween  steel  balls  deviates  greatly  from  the  theoretical  prediction  and  depends  strongly 
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on  the  normal  load. 

Although  the  present  numerical  algorithm  can  simulate  sphere  assemblages 
consisting  of  mviltiple- sized  particles,  an  extension  to  include  ellipsoid  particles  is 
called  for  in  order  to  study  the  particle  shape  effects  on  the  mechanical  behavior  as 
well  as  scalar  transport  properties  of  granular  £issemblages.  To  examine  the  Hertzian- 
contact  resistance  theory,  one  should  find  new  type  of  particles  with  better  surface 
smoothness  and  no  surface  resistance  film.  One  of  the  promising  materials  is  ion- 
exchange  beads,  provided  that  they  are  not  too  soft. 


Appendix  B 
Computer  Code 


program  UIIV3D 

c - Variable  definition: 

c  dm — number  of  dimensions 
c  npt — number  of  particles 

c  idim — number  of  equations  in  quasi-linear  system  of  equation 

c  nl — number  of  microcells  in  X.  Y  and  Z  directions 

c  mznc — maximum  number  of  microcells 

c  mxc — maximum  number  of  contacts  a  particle  could  possibly  have 

c  nnb — the  number  of  adjecent  microcells  around  a  microcell 

c  r--radii  of  particles 

c  x,y,z — position  coordinates  of  particles  at  prior  step 

e  xm,ym,zm — position  coordinates  of  particles  given  by  mean  deformation 

c  xt,yt,zt — position  coordinates  of  particles  given  by  "mean  fluctuation" 

c  ux,uy,uz,«x,vy,uz — fluctuation  displacement 

e  uxm,uym,uzm,vxm,uym,vza — mean  displacement 

c  ttxt , uyt , uzt , vxt , uyt , uzt — total  displacement 

c  a — grand  stiffness  matrix 

c  b — vector  of  unbalanced  force 

c  ak,bk — temparory  6  by  6  matrix 

c  xn — correction  to  the  fluctuating  displacement 

c  dxn — increment  of  xn 

c  res — residual  of  the  system  of  quasi-linear  equations 

c  fnO — normal  force  between  particles  at  prior  step 

c  fnl — normal  force  between  particles  at  curreu'  step 

c  ftO — tangential  force  between  particles  at  prior  step 

c  ftl — tangential  xorce  between  particles  at  current  step 

c  tft — total  tangential  force  between  particles  at  current  step 

c  dft — increment  of  tangential  force  between  particles  at  current  step 

c  tfor — total  force  between  particles  at  current  step 

c  sfx,sfy,8fz — total  force  on  each  particle 

c  smx.smy.smz — total  moment  on  each  particle 

c  nij — fabric  tensor 

c  dl — distance  between  adjacent  particles 

c  nx,ny,nz — contact  normal 
c  s — stress  tensor 

c  stmO.stm — deviatoric  strain  rate  tensor 


76 


76 


c  astrn.psatrn — total  atrain 

c  ahaar — ahear  atrain  incraaent 

c  cond — condnctanca  aatriz 

c  brz,bry,brz — nat  unbalancad  currant  dua  to  naan  potantial  gradiant 

c  cur — intarpartlcla  currant 

c  voltl — fluctuating  potantial 

c  c.alf — alfactiva  conductivity  tanaor 

c  ih — nuabar  of  contacta  for  a  particla 

c  adjc — adjacancy  aatriz  of  aicrocall 

c  kzz.kyy.kzz — for  tha  periodicity  of  aiaulation  coll 

c  aapO — aicrocalla  occupied  by  npt  particlaa  at  prior  atap 

c  aap — aicrocalla  occupied  by  npt  particlaa  at  currant  atop 

c  aappartO — particlaa  occupying  tha  aicrocalla  at  prior  atap 

c  aappart — particlaa  occupying  tha  aicrocalla  at  current  atap 

c  kphi.kthata — contact  noraal  diatribution 

c  ctpx , ctpy , ctpz — coordinatea  of  contact  pointa 

c  thkn.thkt — acalad  aagnituda  of  noraal  and  tangential  force 

c  ashO.aahl — randoa  aaquancaa  of  aicrocalla 

c  khlfl,khlf2 — firat  and  aacond  half  of  a  randoa  aequanca 

c - Definition  of  input  variablaa 

c  iprapk — job  option(0  for  packing,  1  for  repacking,  2  for  relaxing, 

c  3  for  deforaation) 

c  iraax — nuabar  of  real izat ions 

c  iatp — nuabar  of  deforaation  steps 

c  outaax — aaxiaua  nuabar  of  outer  loops 

c  inaax~aaxiBUB  nuabar  of  inner  loops 
c  nracordl,nracord2 — nuabar  of  data  points  to  be  saved 

c  shaara — shear  strain  incraaent 

c  sigaaO — controlling  pressure 

c  epsa — a  snail  nuabar 

c  apso — allowance  for  tha  pressure  balance 

c  apsi — allowance  for  tha  force  balance 

c  dan.i — initial  density  of  loose  packing 

c  danO — desired  packing  density 

c  isead — seed  for  randoa  nuabar  generator 

c  stmO — deviatoric  strain  rate  tensor 

c  rl,r2,r3 — particle  radii 

c  ckn,ckt — noraal  and  tangential  contact  stiffness 

c  cs — interparticle  friction 

c  orf — overralaxation  factor 

c  iral — aaxiaua  nuabar  of  iterations  in  relaxation  aethod 

c  nshuf — nuaber  of  shuffling 

c  Iflip — nuaber  of  flip  in  shuffling 

c  nrif — nuaber  of  riffle  in  shuffling 

c  itrplaax — aaxiaua  nuaber  of  trials  in  placing  a  particle  in  a  aicrocell 

c  without  ovarlaping  adjacent  particles 

c  daval  to  dv8 — controling  paraaaters  for  packing 
c  i_cuau — an  integer  nuabar 

c  itriala — aaxiaua  nuabar  of  trials  of  packing  to  desired  density 

c  a_ya — elastic  aodulus 


pois_r — Poisson'i  ratio 

barkn — acaliiig  factor  of  noraal  contact  atiffnass 
barR — scaling  factor  for  radii 
grad — potential  gradient 

conk,conb — two  paraneters  in  load~resistance  relation 
isign — pressure  control  paraaater 


iaplicit  real  (a-h,o-z) 
iaplicit  integer  (i-n) 
integer  dm 

paraaieter(dBs3,npts48,idiB=6*npt,nl=8.amc=nl*edB,Bucc=12) 
paraaeter  (dns3  ,npt=48,  idisFdenpt  .nlslO  ,aucnc=nle*dB,Bucc=12) 
paraaeter  (dBs3 ,  npt=100 .  idia=6*npt  ,nls9  ,axnc=nl*eda,azcsl2) 
paraaet  er (nnbs342 , pi=3 . 1415927 , zero=0 . 0 , izero^O , ones 1 . 0 , Tal= 1000 . ) 
real  x(npt) ,y(npt) ,z(npt) .r(npt) 

real  xa(npt) ,yB(npt) ,za(npt) ,xt(npt) .yt(npt) ,zt(npt) 
real  ux(npt) ,ay (npt) ,nz(npt) .uxaCnpt} ,uya(npt) ,nza(npt) 
real  Bx(npt) ,«y (npt) ,vz(npt) .«xa(npt) ,sya(npt) .«za(npt) 
real  nxt(npt) ,uyt(npt) ,azt(npt),nt(npt).syt(npt) ,vzt(npt) 
real  a(idia,idin) ,b(idin) ,ak(6,6),bk(6,6) ,xn(idia) 

real  fnO(npt ,npt) ,fnl(npt,npt) ,tft(npt,npt) 

real  ftO(npt,npt,da) ,ftl(npt,npt,da) ,dft(da) 

real  sfx(npt) ,sfy(npt) ,sfz(npt) ,sax(npt) ,say(npt) ,saz(npt) 

real  tf or(npt ,npt) ,dl(npt ,npt) ,dxn(idia) 

real  nx (npt , npt ) , ny (npt , npt ) , nz(npt , npt) , ni j (da , da) 

real  ras(idia) ,8(da,da) 

real  ctpx(npt,axc) ,ctpy(npt,BXc) ,ctpz(npt,nxc) 

real  8tm0(3,3) ,8tm(3,3) ,8stm(3,3),p8Stm(3,3) ,8tarin(3,3) 

real  cond(npt ,npt) ,brx(npt) ,bry(npt) ,brz(npt) 
real  cur(npt,npt) ,voltf (npt),c_eff (3,3) 
real  cxn(npt) ,cre8(npt) ,cdzn(npt) 

integer  ih (npt ) , kxx (axnc , nnb) , kyy (axnc , nnb) , kzz (axnc , nnb) 

integer  ad jc (axnc, nnb) ,Bcont(npt,npt) ,Bap(npt) ,mapO(npt) 

integer  Bappart(nxnc) ,aappartO(Bxnc) 

integer  8test,outBax,high8tm,kphi(12) ,ktheta(24) 

integer  thkn(npt,Bzc),thkt(npt,Bxc) 

integer  BshO(axnc)  ,ashl(Bxnc)  ,khlf  i(Bznc)  ,Uilf2(BUcnc) 

open (units 10 , f iie= • 3dat ' , statuss  *  old ’ ) 

open  (unit =30 , f iles ' condu  * , 8tatu8= ’ unknom ' ) 

open (unit =49 , f ile= ’ inner ’ , 8tatU8= ' unknown ' ) 
open (unit =50 , f ile= ’ outer  * , status= ’ unknown ' ) 
open (unit=5 1 , f ile= ’ displ ‘ , 8tatus= 'unknown* ) 
open(unit=53 , f ile= ' slope ' ,stattt8= ' unknown ' ) 
open (unit=54 , f ile= ' phi ' , 8tatus= 'unknown ' ) 
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op«n(iuiits55 ,  i  il**  *  nbond  * ,  stata«= '  nnlrnom ' ) 
op«ii(uBitsS6,f  lies*  theta*  .statna^’anknoni') 

op«n(unit«61 ,lil«»*pT«datO* .statvsB’nnknoni' ) 
open  (unit  s62 , 1  ilas  *  coordO  * ,  status  *  naknoni  * ) 
opsxi(iuiits7 1 ,  f  ilss  *  prsdat  1  * ,  atatass  *  nnknovn  * ) 
open  (units72 , 1  lias '  coord  1  * ,  status '  uknovn  * ) 

opan(uit368 ,  i ilas  *  luvO ' ,  status  'unknown  * ) 
opan(uits69 ,  files  >  f  nav  * ,  status*  unknown' ) 
open (uits70  ,f  lies  *  f  n* ,  statuss  'unknown ' ) 

open (unit s80 .files \coord ' .statuss 'unknown ' ) 
open(unit=81 .files  >  contc  * . statuss 'unknown ' ) 
open  (uit  s82 .  f  iles  *  contp  * .  statuss '  unknown  * ) 
open  (unit =83  .files  •  thick ' .  status  'unknown ' ) 
open (unit  s92 . f iles  *  fore  e ' . statuss ' unknown  * ) 
c  open(unit=98.f iles'balan' .statuss 'unknown') 
open (unit =99 . f iles ' probl esi  * . statuss ' unknown ' ) 

c - reading  initial  data  and  parameters 

read( 10 .*) iprepk. irmax. istp.outmax. inmaz .nrecordl .nrecord2 
readdO.*)  shearm .  a  igmaO .  epsa .  epso .  epsi .  den_i .  denO .  is  e  ed 
read(10.*)8tm0(l.l)  .8tm0(l,2).8tm0(1.3) 
read(10.«)stm0(2.1)  .8tm0(2.2).stm0(2.3) 
read(10.e)8tm0(3.1)  .8tm0(3.2).stm0(3.3) 
read(10.*)rl.r2.r3.ckn.ckt.cs.orf .irel 

r ead ( 10 . * } nshuf . If lip , nr if . itrplmax . devnl . devm2 . devB3 . af 1 . af 2 

r  ead  ( 1 0  .<•>)  den02 .  denO  1 .  deni .  den2 .  den3 .  den4 .  dens 

read(10.«)dv0.dvl.dw2.dv3.dv4.dv&.dv6.dv7.dv8 

read( 10 . e) i.eumu. itr iala . e_ym.pois_r.barkn.barR , grad . amda 

read( 10 . e) conk . conb . isigm 

write(e.e} 'isigm' .isigm 

if (nrecordl . ne . 0}then 
idistlsistp/nrecordl 

else 

idistlslOOOO 
end  if 

if (nrecord2 . ne . 0}then 
idist2=i8tp/nrecord2 

else 

idi8t2s 10000 
end  if 

if (istp . le .nrecordl) idist 1=1 
if (istp . le . nrecord2) idist2=l 
monops 1 

if (rl . ne . r2 . or . r2 . ne . r3)monop=2 


if (ipr epk . le . 2) then 


•pai>0.0 

cssQ.O 

call  rinit2(3.3,strn0.0.0) 
if  (danO .  la .  0 . 6)  shaansQ .  0 
and  if 

- carrying  out  "imax"  raalizations 

do  599  irsl.iraaz 
itrial^l 
199  icnausO 

- ganarating  initial  loosa  systan  for  packing 

if ( ipr apk . aq . 0) than 

iatoldsQ 

avsO.O 

n203nint  (raal(npt)  /  ( 1 . 0+(r2/rl )  aadn-t-  (r2/r3)  **da)  ) 

nlOsnint (raal(n20) * ( (r2/rl }aadB) ) 

n30=npt-nl0-n20 

n:ita(*,*)nl0.n20,n30 

pl=raal(nlO)/raal(npt) 

p2sraal  (nl0-«-n20)  /r  aal  (npt) 

solidvs4.0api*(nl0*rl**dJB'i-n20*r2**dm'4-n30ar3**dB)/3.0 

i  s  aads  (is  aad-»’int  ( 1 00000*ran(  i8ead))}*2-fl 

ratioxz=1.0 

ratioyz=1.0 

stmaxsahaara 

dalz0=2 . 0*rl/8qrt ( (ratioxz+tanCatmax) ) ♦♦2+ratioyz**2+l . 0) 
dalxO=dalzOarat ioxz 
d8ly0=dalz0*ratioyz 
nn=l'*-axp(log(8olidT/dan_i)/da)/d8lz0 
101  vrita(*, a) 'initial  na: '  ,iu> 

103  if  (na.gt.nDthan 

writaC*,*) ’nB>nl,not  anongh  aicrocall. 'nn,nl=' , na.nl 
arita(99,a} 'nB>nl,aot  anough  aicrocall. 'na.nls' .na.nl 
atop 
and  if 

dxO=dalxOann 

dyO=dalyOana 

dzO^dalzOana 

ncall^nnaada 

i  j  80=l-i-2 .  Oar3/axp(log(8olidv/(ncall*denO)  )/dB) 
if (raal (i j  sO) . ga . raal (na)/2 . 0) then 
nasnaal 
go  to  101 


80 


•nd  it 

call  adj3(aznc.lja0.iui.iuib,adjc,kxx,ky7.kzz) 
call  Bhuttla(ncall,nshiit,ltlip,iiLrif ,isaad,uhO. 

aahl.khltl, kbits) 

call  vac_$iinit(BappareO,Bxac,izaro) 

al^O 

n2-0 

n3-0 

kapsQ 

BcountsO 

do  100  Bpsl.npt 

it (aonop. aq. l)tb«n 
r(np)=r2 

alaa 

arita(*,a) 'pleas#  include  T_aasign  subroutine' 

call  r_assign(nl.n2.n3,nl0.n20.n30,rl.r2.r3, 

+  pl,p2.iaeed.r(np)) 

end  it 

83  kap=kap'fl 

it (kap . gt . ncell ) then 

c  urite(*,*)'ir=',ir,'  n«»',n«,'  np=',np,'  kap=',kap 

nnsxm-t'l 
go  to  103 
end  it 

klsnshO(kap) 

ka=kl/(nmenn) 

it  (kBeniBenB.ne.kl)kB=]Qa'i'l 

jB= (kl- (km-1) ennenm) /na 

it  (  jaena.na .  (kl-(kB-l)*naenn)  )  ja=  ja-i-1 

iaskl- (ka-l ) enaena-C ja-1) sna 

ncount=ncount-i- 1 

hx= ( ia- 1 ) edelzO 
hy=(ja-l)*delyO 
hza (ka- 1 } sdelzO 

itrpl*0 

84  itrpl=itrpl+l 
z(np}=hz-t’ran(i8eed)edelx0 
y  (np)  ::hyi-ran(iseed)  edelyO 
z  (np)  =hz-i-ran(  iseed)  edelzO 

do  80  j=l,nnb 

kn^aappartO (adj  c(kl , j ) ) 


il(kB.*q.O)go  to  80 
k2>»ap0(lai) 

zknsz (ka) -kxx (k 1 . j ) *dx0 
7kn*j(kB)-k77(kl , j)*d70 
zka»z (kn) -kzz (kl . j ) *dz0 

d02« (z (up) -zka) ♦ *2+ (7 (np) -pka) **2+ (z (ap) -zka) **2 
il(d02.go.(r(np)-»^r(ka))*(r(np)'«'r(ka)))go  to  80 
if  (Itrpl .  It .  itrplBaz)tliea 
go  to  84 

olao 

go  to  83 
and  if 

80  continaa 

■ap0(np)skl 

BappartO(kl)=np 

100  continna 


call  rinit2(npt,npt,fn0,zero) 
call  riait3(apt,Bpt,dB,ftO,zaro) 
c  «Tita(a,a) •acouat** ,acoaat 

daaaitpssolidv/CdzOadpOadzO) 

vrit a (62 , a) ir , iat , dzO .dpO ,dzO , aa, i j sO 
do  79  i=l,apt 

writa(62,*)i.r(i),z(i),7(i),z(i) 

79  coatiaua 

c - atartiag  froa  a  packad  87ataa  to  rapack,  to  ralaz  tha  atraaa, 

c  to  do  aiapla  ahaar  ate. 

alaa 

vrita(*,*)'atart  froa  packad  a7ataa* 
call  riait2(npt,apt,fa0,zaro) 
call  riait3(npt,npt,da,ft0,zaro) 

raad(61 ,a)irt,nl0,a20,o30,pl,p2,dan8it7,aolidv,8igBa 
raad(61,a)ia,i8told,88tra(l,l) ,88tra(l,2)  ,88tm(l,3) , 
a8tra(2,l) ,88tra(2,2),a8tza(2,3) , 
a8tra(3 , 1 ) , 88tra(3 , 2) , 88tzn(3 , 3) , av ,dalav_la8t , ia idaO 

do  149  i=l,apt 

raad(61 , *) it .ataa 
if (atoa.oq.O)go  to  149 
do  147  j=l,atoa 

raad(61.*)jla.fo0(i,jla).ft0(i. jla.l), 
ft0(i.jla,2),ft0(i,jla,3) 
coatiaua 
coatiaua 


+ 

+ 


147 

149 
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r«ad(62 ,«) irt , istt ,dxO ,dyO .dzO .amt . i jsOt 
do  176  i>l,apt 

r«ad(62.*)it,r(i),x(i).j(i),z(i) 

contiauo 

if ( ipr opk . oq . 1 ) then 
denspsdeniity 
deidea’’ (den0-denap)/10. 0 
nden-1 
end  if 

ratioxz-dxO/dzO 
ratioyzsdyO/dzO 
if (iprepk . eq. 3)then 

etnaxsnax  (  aba  (  stxnO  ( 1 , 2)  ) ,  abe  (  etmO  ( 1 , 3  )  ) ,  abs  (  8  trnO  (  2 , 3)  ) , 
abs(stm0(2,l)),aba(strn0(3.1)).aba(8trn0(3.2)))*(istold-t-istp) 

else 

atznax-O . 4 
end  if 

delz0s2.0erl/8qrt((ratioxZ'«'tan(8tZBax))**2-fratioyzee2-t'l  .0) 
nn=dzO/delzO-*'  1 

if (na . gt . nl) then 

vrite(e,*) 'n>>nl,not  enough  nicrocell. ’ , 'aa.nls’ .na.nl 
eTite(99,«) 'nn>nl,not  enough  nicrocell. *nm,nl=' , nn.nl 
atop 
end  if 

oriteC*.*)* computed  tut*. an,*  nl'.nl 
atop 

na=nl 

delzO=dzO/nn 

delxO=delzOeratioxz 

delyO^delzOeratioyz 

ncell=naeedn 

i  j80’:l-*’2 .  Oer3/exp(log(8olidv/(ncelledenO)  }/dn) 
if (real(ij80).ge.real(nn)/2.0)then 
nn=nn-t-l 
go  to  177 
end  if 

write(*,*)*nn«*,nn,»,  ijaOs* .ijaO, * ,  ncell=* .ncell 

call  adj3(nxnc,ij80,nn,nnb,adjc,kxx,kyy,kzz) 
call  napping (apt , nn ,nxnc , aatzn , delxO , delyO , delzO , 

X , y , z , napO , nappart  0 , iout a ) 

end  if 

kf lag-0 
VTite(98,*)ir 


call  riait2(3,3,s>trB.0.0) 
call  riBit2(3,3,ptstrn,0.0) 


c 

c — a 


692 


690 


691 

693 


ilargsl 

iiaalsO 

iX ( ipr apk . aq.  2) istps 1 

number  of  cyclic  sbaaringa  for  packing  or  stapa  for  sinpla  shaar 
do  598  istsistold-t'l.istold-fistp 
icback=0 


itar3s0 

loustmsO 

highatm-O 

noatop^O 

if ( itar3 . no . 3 . and . noatop. aq. 0)than 
ichack^ichack^l 

call  niatx_copy(3,3,etm0,atm) 

alaa 

if  (loaatm .  It .  2  )  than 
ichacksichack-t'l 
do  690  i-1.3 
do  690  j=l,3 
8tm(i,  j)=8tm(i,  j)/2.0 
continue 

lovatmslouatm-M 

itar3=0 

if (noatop . aq. I)no8top=0 

alaa 

if  (highatm .  aq.  0)then 

ichack^ichack-fl 
do  691  i=l,3 

do  691  j=l,3 

8tm(i,  j)=8tm0(i,j)*2.0 
continue 

alaa 

icheck=ichack+l 
do  693  i=l,3 

do  693  j=l,3 

8tm(i,  j)=8tm(i,j)*2.0 
continue 

end  if 


if (noatop . aq . 1 ) no8top=0 

itar3=0 


high8tm=high8tm-i'l 
if  (highatm .  gt .  2 )than 


«rita(99,«) 'aorry,  I  caa&ot  find  a  satisfactory  dslsv' 
writs (99,*) 'to  balaacs  tbs  strsss,  stop!  highstm' 
stop 

snd  if 
snd  if 
and  if 

if  (iprspk .  Is .  1  )thaii 

if (dsnO . Is . 0 . 6 . or . (dsnO . gt . 0 . 6 . and . dans ity . gt . danO-O . 002 
.  and .  abs  (psstm  (3 , 1 )  -t^psatm  ( 1 , 2)  ) .  It .  1 .  Os-8)  )  than 
stm(3,l)s0.0 
stm(l,2)*0.0 

slss 

ismal^isnal-fl 
if (■od(ilarg,2) .ns.0)thsn 

if  (ismal.  sq.  l)8tm(3,  l)sshsarB 
if  (isaal .  aq .  2)  stm(3 , 1 )  =-shaazB 
if (i8nal.aq.3)stin(3,l)*-8bsaxB 
if  (isnal .  aq  .4)  stm(3 . 1)  *8hsarm 
8tm(l,2)=0.0 
88tm(l,2)*0.0 
slss 

if  (isnal .  aq .  1 ) stm(  1 , 2) =sbsam 
if  (  isnal .  aq .  2)  stm(  1 , 2)  *-shaam 
if  (isnal .  aq .  3)  8tm(  1 , 2)  *-8haam 
if  (isnal .  aq . 4)  stm(  1 , 2)  *shaam 
stm(3,l)=0.0 
sstm(3,l)=0.0 
and  if 

if (nod ( isnal , 4) . aq. 0) than 
ilarg^ilarg-t-l 
isnal^O 
and  if 

and  if 
snd  if 

if (ichack .gt . 12)than 

writa(99,*) 'progran  teminatsd  das  to  unconvarged  outloop' 
writa(6,*) 'progran  taminatsd  das  to  onconvargsd  oatloop* 

stop 
and  if 

8tm_c=0.0 
88tm_c*0.0 
do  687  i*l,3 
do  687  j=1.3 

88tm(i,  j)*p88tm(i,  j)'fstm(i,  j) 


il (&bs (atraCi . j ) ) • gt . abs (stni_c) ) stni_csstrn( i , j ) 
if  (abs (■•trad ,  j  ))  .gt  .abB(s8tra_c) )  sstra.cssstrad ,  j  ) 

687  coatiaus 

c  vrits(99,s)ir.ist,it«Toat,stTa_c 

iter out sO 

c  vrite(SO ,*)‘ ir , iat . iterout . stest .test .error .del ev . ev .f In/a. 

c  -t’deas  * 

«rite(e.*)*ir.i8t. iterout .stest. test. error . dele v . ev . i lu/m , 
-t'deas.  isiga' 

uriteC*,*)'  • 

c  urite(98.*)*  * 


c - choosiag  the  iacreaeat  ol  coatractioa  lor  packiag 

il (iprei^. le . l)thea 

it (iprepk . eq. 0 . aad. ist . le . l)thea 
delevO*-dvO 

else  il(dea8ity.ge.dea0)thea 
delevQsO . 0 

iCUBU^iCUBU-t-l 

else 

il (deas itj . It . dea02) thea 
delevO=-dvl 

else  il (deasity .It .deaOl . aad. deasity . ge . dan02) thea 
delev0»-dv2 

else  il (deasity . It . deal . aad .deas ity . ge . deaOl ) thea 
delev0=-dv3 

else  il (deas ity . It . dea2 . aad .deas ity . ge . deal ) thea 
delev0=-dv4 

else  iKdeasity. It. deaS. aad. deasity. ge.dea2)thea 
delev0=-dv5 

else  il (deasity .It. dea4. aad. deasity .ge.deaS) thea 
delev0=-dv6 

else  il (deas ity . It . deaS . aad . deasity . ge . dea4) thea 
delev0=-dv7 
else 

delev0=-dv8 
ead  il 

ead  il 

delev=delevO 

c - do  packiag.  lorce  balaaciag,  aad  computiag  stress 

call  delev.stress (da, noaop. apt .idia. am. iprepk. ir. ist .deaO. 

+  iterout , iaaer , iter2 , it er3 . idistl , zero , val , cka , ckt , cs , 

+  iaazLX , epsi , epso , epsa, irel , sigaaO , isiga.laaverage .laav.cur . 

*  delev , ev , sstra , stra , delxO .delyO , delzO , delz , dely , delz , axac , aab , 
-*■  r.x.y.z.za.ya.za.xt.yt.zt.uz.uy.uz.uza.uya.uza.uxt.uyt.nzt, 

+  ux,uy,uz,«xa,«yB,uza,uxt,uyt,0zt,8olidv,voluae, deasity, 

+  laO , lai , ItO ,lt 1 , dlt ,dl , six , sly , slz , sax , say , saz , itrial , 
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*  a,b,ak,bk,adjc,kxz.k77>)^i>^>&Jin2.MpO,Kap.nbon.kslp,kact. 
za,dzn,r«s,orf .nuBT.BappartO.aappart. 

*  a .  t  aat ,  s  igna .  staat ,  error .  errontia .  delevain .  kllag ,  aada) 

go  to  604 
end  11 

c - algoritlue  to  achieve  the  desired  isotropic  stress 

iter2=l 

iter3=0 

inega=l 

i¥ant=0 


devl=zero 
delev=devl 
if ( ipr epk . eq . 2 ) then 
c  iteroutaiterottt-t^l 

c  test=sigBa-8ignaO 

c  l_devl=test 

c  il(te8t.lt.0.0)then 

c  8teat=-l 

c  else 

c  stestsl 

c  end  il 

c  error^test/signaO 

else 

iterout«iterout+ 1 

call  delev.stress (dn , nonop , npt , idin ,  nn , ipr epk , ir , ist , denO , 
iterout ,  inner ,  iter2 ,  iter3 ,  idistl  .zero .  val ,  ckn ,  ckt ,  cs , 
innu ,  eps i ,  epso , epsa , ir el , s ignaO ,  is igm .Inaverage , Inav.cur , 

+  delev . ev , sstm , strn , delzO .delyO . delzO , delz , dely , delz , mznc , nnb , 
r.z.y.z.zn.yn.zm.zt.yt.zt.uz.ay.uz.uzn.uyn.uza.uzt.uyt.uzt. 

+  wz , vy , wz . wzB , vym , szn , szt . vy t . vzt . sol idv , volume .density, 

*  InO.lnl.ltO.ltl.dlt.dl.slz.sly.slz.smz.smy.snz.itrial, 

a,b,ak,bk,adjc,kzz,kyy.kzz,nz,ny,nz,mapO,map,nbon,kslp,kact, 

-t-  zn,dzn,res,orl,numr,mappart0.mappart. 

+  8 ,  t  est ,  s  igna ,  st  est ,  error ,  erromin ,  delevnin ,  kl  lag ,  anda) 
il(abs ( error ). It. epso)go  to  604 

l_devl=te8t 
end  il 
iside=8te8t 
testO^test 
errorO=error 

S92  igo942=0 

il(isideeinega.eq.-l)then 

il (ist . le . 1 . or . ivent . eq. l)then 
devn_l=devnl 
else  il(iter3.ge. l)then 

devn_l=delevmin-0. 0005*iter3 


•It* 

if (iaidaO . •q.O)th«n 
davm_lsdaTBi 

elsa  if (iBidaO.aq.-l)tliaii 
da VB_ 1 saf l*dalaT_laat 

alia 

davB_ls-af l*dalaT_laat*4 . 0 
and  if 
and  if 

dalav=daTB_l 

itaroutsitarout'fl 

call  dalav.atrasa (da , aonop , npt , idia, na , ipr apk , ir , ist , danO , 

*  itaront , innar , itar2 , itar3 , idiatl , zaro , val , ckn, ckt , ca . 

+  inaaz , apai , apao , apaa, iral.aigaaO , iaiga.f navaraga ,f nav.cnr , 
dala  V ,  av ,  aatrn .  atrn ,  dalzO  .dalyO ,  dalzO ,  dalz ,  daly ,  dalz ,  axnc ,  nnb 
r,z,y,z,za.yB,za,zt,yt,zt.uz,Qy,nz,nza,nya,uza,iizt,uyt,uzt, 

*  vz.vy ,vz,na,vya,vza,«zt.vyt,vzt,aolidT,voliiaa,danBity, 
f  nO ,  f  nl ,  f  to  ,f  1 1  ,df  t  ,dl ,  af  z .  af  y .  af  z ,  aaz ,  aay ,  aaz ,  itrial , 

+  a,b,ak,bk,adjc,kzz,kyy,kzz,nz,ny,nz,aap0,aap,nbon,k8lp,kact, 

*  zn,dzn,raa,orf ,nuar,BappartO,aappart, 

a ,  taat ,  ai^pna,  ateat ,  azror ,  arrozain  .dalavain ,  kf  lag ,  aaula) 
if (abaCarror) .lt.apao)go  to  604 

favmlateat 

if (iat . la . 1 . or . ivant . aq. l)tkan 
davB_2sdavB2 

ivantsO 

alaa  if (itarS.ga.Dthan 

davBi_2=dalevBdtt~0.000laitar3 

alaa 

if ( ia idaO . aq . 0) than 
devB_2=davn2 

alaa  if (iaidaO.aq.-Dthan 
davB_2=af2adalav_la8t 
alaa 

davB_2=-af 2*dalav_laat*4 . 0 
and  if 
and  if 

dalav=davn_2 
it  arontsit  arout-*- 1 

call  dalav.atraaa (dm , monop , npt , idim , nm , ipr apk . ir , ia t , danO , 

-I-  iterout ,  innar ,  it ar2 ,  it ar3 ,  idiat  1 ,  zero ,  val ,  ckn ,  ckt ,  ca , 

-i-  inmaz ,  apai , apao , apaa, iral ,aigmaO , iaigm ,f naver2ige ,f nav.cur , 
delav ,  av ,  aatm ,  atm , dalzO, dalyO, dalzO , dalz , daly , dalz , mznc , nnb 
-*■  r,z,y,z,zm,ym,zm,zt,yt,zt,nz,iiy,nz,uzm,nym,nzm,uzt,nyt,uzt , 
vz,vy,vz,vzm,vym,vzm,vzt,vyt,vzt,8olidv, volume, denaity, 

*  fnO,fnl,ftO,ftl,dft,dl,8fz,afy,afz,amz,amy,amz, itrial, 
a,b,ak,bk,adjc,kzz,kyy,kzz,nz,ny,nz,mapO,map,nbon,kalp,kact, 
zn,dzn,raa,orf  ,nnmr,mappartO,mappart, 

a ,  taat ,  a  ipaa ,  ataat ,  error ,  arrormin ,  dalavmin ,  kf  lag ,  amda) 


if (abs(«rTor).lt.«pao)go  to  604 


f«ni2-t«st 
go  to  791 

•ls« 

942  if ( ( (iaidoO .!•. 0 .or. arrorO . It . 1 . 0) . and. iterS . aq.O) . 

or.igo942.eq.l)than 

davB_2-0.0 

favB23taatO 

deni_lsdavB3 

dalav=davai_l 

941  itaroutsiterout't'l 

call  delav_atra8a(dn,Bonop,npt,idiB,iuB,iprepk,ir,iat,denO, 
itaront ,  inaar ,  itar2 ,  it  ar3 ,  idiatl ,  zaro ,  val ,  ckn ,  ckt ,  ca , 
inmax ,  apa  i ,  apao ,  apaa,  iral ,  aignaO .  iaign ,  f  navaraga ,  f  nav_cur , 

*  dalav ,  at ,  aa tm ,  atm ,  dalxO .  dalyO .  dalzO ,  dalx ,  daly ,  dalz ,  mxnc ,  nnb , 

*  r,x,y,z,xm,yn,zm,xt,yt,zt.tix,ny.az,nxm,nym,uzm,uxt,uyt,uzt, 

+  vx,vy,vz,vxn,vym,vzn,Bxt,vyt,azt,aolidv,voluma,dan8ity, 

fn0,fnl,ft0,ftl,dft.dl,8fx,8fy,afz,  anx ,  any ,  aaz ,  itr  ial , 
a,b,ak,bk,adjc,kxx,kyy,kzz,nx,ny,nz,napO,map,nbon,kalp,kact, 
xn,dxn,raa,orf  ,nnnr,nappartO,nappart, 

-i-  a ,  t aa t ,  a igma ,  at aat ,  arror ,  arrormin , dalavmin ,  kf lag ,  amda) 
if (ab8(arror).lt.apao)go  to  604 

f avml^taat 

if (favnl.gt.favB2)than 
go  to  791 
alaa 

davm_l=10.0*davm_l 

dalav=d8vm_l 

go  to  941 
and  if 
alaa 


if (it  er3 .aq.O) than 

davn_l=dalav_la8t/af 1 
alaa 

da vm_ 1 =dala vnin/ (afiaitar3) 
and  if 

dalav=davn_l 

itarout^itaront't’l 

call  dalav.atraaa (dm , nonop , npt , idin , nn, ipr apk , ir , iat , danO , 
itaront ,  innar ,  itar2 ,  itar3,  idiatl  ,zaro ,  val ,  ckn,  ckt ,  ca , 

+  innax , apai , apao , apaa, iral , a ipaaO, iaign ,f navaraga .fnav.cur , 

+  dalav ,  av ,  aatm ,  atm ,  dalxO  ,daly0 ,  dalzO ,  dalx ,  daly ,  dalz ,  nxnc ,  nnb , 
+  r,x.y,z,xn,yn,zm,xt,yt,zt,nx,ny,nz,nxn,nyn,nza,uxt,uyt,uzt, 
«x,vy,vz,«xn,vyn,vzn,vxt,vyt,«zt,8olidv,volnna,danaity, 
■•■fn0,fnl,ft0,ftl,dft,dl,8fx,afy,  af  z ,  anx ,  any ,  anz ,  itr  ial , 

a,b,ak,bk,adjc,kxx,kyy,kzz,nx,ny,nz,nap0,nap,nbon,k8lp,kact, 

*  xn,dm,ra8,orf ,nunr,nappartO,nappart, 
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+  s ,  t«st  ,»igui,  «tMt ,  urror ,  •rrond.a  .dalcvBin ,  kllag ,  awla) 
il(abs(«rror).lt.*pso)go  to  604 

lOVBlstOtt 

il(it«r3.oq.0)thoB 

dovB_2>d6lov_last/al2 

•Iso 

d«VB_23d«l«TBia/ (s22* it orS) 

•nd  ii 

d*l«vsd«va_2 

it«routsit«rottt-*-l 

call  del«v_str«aa(dB.Bonop.npt , idia.ui,  iprapk,  ir ,  ist  ,daa0, 

*  itorout ,  iim«r ,  lt«r2 ,  it  arS ,  idist  1 ,  zaro ,  val ,  ckn ,  ckt ,  cs , 

+  laaiaz , apai , opao , apaa, ir al , aigaaO . iaija , Inavoraga , f nav.cnr , 
dalov ,  OT ,  aatm .  at  m ,  dolzO  .dolyO ,  dolzO ,  dolz ,  daly ,  dolz ,  aocnc ,  nnb , 

*  r,z,y,z,ZB,ya,za,xt,yt,zt,ttz,ay,iiz.uai,nyB,iizB.axt,iiyt.uzt, 

+  vz,«y,«z,nai,«ya,«ZB,«zt,«yt,azt,aolidv,voliaia,donaity, 

InO ,  f  nl ,  ItO  ,ttl  ,dlt  ,dl .  alz,  af  y  ,alz ,  aaz ,  aay ,  aaz ,  itrial , 
a,b,ak,bk,adje,kzz,kyy.kzz,nz,By,nz,aapO,Bap,nboii,kslp,kact, 
zn.dzn.raa.orl.niwr.BappartO.aappart, 
a ,  t  aat ,  a  ipaa ,  at  aat ,  azror ,  orrozain  ,dola  vain ,  kf  lag .  aada) 
il(abs(orror).lt.apso)go  to  604 

fova2staat 

il (t aval . gt . fova2. and . f aval . gt . 0 . 0)than 
go  to  791 

also 

ii (taatO . gt . f •va2} than 
iavalstoatO 
dava_l~0.0 
go  to  791 

also  if (taatO.gt.fovaDthon 
fava2-faval 
dava_2=dava_l 
f aval=tastO 
dava_l=0.0 
go  to  791 
•Isa 

igo942=l 

go  to  942 
and  if 
and  if 


•nd  if 
and  if 

791  if (favBl*fava2.1t.0.0)than 

if (faval.gt.O.O)tban 
daTl»davB_2 


f_deTl=lam2 

else 

devl^dan.! 
d«vh>d«Ta_2 
f_d«vl>favBl 
f_devhsf«vB2 
end  il 
go  to  940 
end  11 

slpks(leral-leYn2)/(devB_l-deTB_2) 
if (abs(slpk) .lt.0.02)then 
if (i8ide«inega.eq.-l)then 
dnvB_lsdevB_l-0 . 002 
delav=devM_l 

it  erout «it  eront+l 

if (it er2 . eq . 1 . and . it eront . gt . 30) then 
nostop=l 
go  to  692 
end  if 

call  delev.stress (da .nonop.npt . idia.na,  iprepk ,  ir ,  ist , denO , 

*  iterout , inner , iter2 , it erS , idiatl , zero , val , ckn , ckt , c> , 

+  inaax , epsi , epso , epaa, irel.eigaaO , iaiga ,f narerage ,f nat.cnr , 

*  delev, ev.aatm, atm, delxO,delyO,delzO,delx»dely,delz,nmetimb, 

*  r,x,y.z,xa,yn,za,xt,yt,zt,uz,ny,nz,nra,nya,uza,ttxt,nyt,nzt, 

+  vx ,  vy ,  vz  ,«xa,  vyn,«za,  vxt  ,«yt  ,«zt ,  aolidv ,  volnae ,  denaity , 

+  fnO,fnl,ftO,ftl,dft,dl,afx,afy,afz, aax , any , aaz , itr ial , 

*  a,b,ak,bk,adjc,kxx,kyy,kzz,nx,ny,nz,BapO,nap,nbon,kalp,kact, 
xn,dxn,rea,orf ,nnar,BappartO,Bappart, 

+  a , t  eat , a igaa , at eat , error . err ozain , delevain , kf lag , aada) 
if (aba(error) .lt.epao)go  to  604 

feval=teat 
go  to  791 
elae 

if (deva_l . le . 0 . 0)then 

deva_l=2 . 0*deva_l 

elae 

deva_l=0. 5edeva_l 

end  if 

delev=devB_l 

iterout=iteroat-t'l 

if ( it er2 . eq . 1 . and . it eront . gt . 30 ) then 
no8top=l 
go  to  692 

end  if 

call  delev_atre8a(dB,Bonop,npt,idiB,nB,ipropk,ir, iat.denO, 
iterout ,  inner ,  iter2 ,  itorS ,  idiatl ,  zero ,  val ,  ckn ,  ckt ,  ca , 
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*  iiuuuc,«psi,«pso,«ps«,ir«l,iipMLO,isi{B,fiiav«raf«,fB&T_cur, 

-t'  d«l«Y ,  CT .  sstru ,  strn .  dalxO ,  d«lyO ,  dalzO ,  d«lz ,  d«ly ,  d«lz .  wmc ,  anb , 

*  r,z.7,z,za,yB,ZB.zt,yt,zt,iiz.ny,iiz,iuB,nyB,iizB,iirt,nyt,uzt, 

*  vz,zy,«z,«XB,vyB,vzB,nt,«yt,nt,zolidT,TolnB«,d«a«ity, 

*  faO,fnl,ltO,ltl.dlt,dl.ilz,sly.slz,uu,My,Baz,itrial, 

*  z,b,zk,bk,adjc,kzz,kyy,kzz,nx,Ay,iiz,upO.Bap,nboB,kslp,kact, 
ZB,dzii.r«z,orl,iuiBr,aappartO,Bi^art, 

*  s ,  t ast .  sigBZ ,  ztMt ,  azTor ,  arrozBia  .dalarain ,  kf lag ,  aada) 

il(aba(arTor).lt.apso)go  to  604 

favBlstast 

go  to  791 
aad  11 

alaa 

ilovsO 

iliigli=0 

davaL3=davB_l-laval/slpk 
il(isida*iBaga.aq.-l)thaa 
11 (davB_3 . It . 0 . Olthaa 
dalaT=dara_3 

also 

dalavs-O.OOOOS 
and  11 

dav»_3=dalav 

alsa 

11 (davn_3 . gt . 0 . 001 )dava_3«0 . 0005 
dalavsdava_3 
and  11 

ltaroats=ltarout+l 

ll(ltar2.aq.l. and . It arout . gt . 30 ) than 
noatop-1 
go  to  692 

and  11 

call  dalav.strasa (da , nonop, npt , Idla.nn, Iprapk ,  Ir ,  lat , danO , 
Itarottt ,  Inner ,  ltar2 ,  ltar3 ,  Idlstl ,  zaro ,  val ,  ckn ,  ckt ,  ca , 

*  Innax , aps 1 , apso , apa  a , Ir al . a IgnaO , la Ign , Inavaraga , lnav_cnr , 

+  dele  V ,  av ,  aatm ,  atm ,  dalzO ,  dalyO  ,dalzO ,  dalz ,  daly ,  dalz ,  nznc ,  nnb , 
+  r.z.y.z.XB,yn,zn,zt.yt,zt,nx.uy,uz,iuai,nya,uzB.nzt,uyt.nzt, 

+  vx,vy,az,vzB,vyB,vza,¥xt,«yt,vzt,aolldv,TolnBa,danalty, 
lnO,lnl,ltO,ltl,dlt,dl,alx,aly,alz,aBX,aBy,anz,ltrlid, 
a,b,ak,bk,adje,kxx,kyy,kzz,nx,ny,nz,BapO,Bap,nbon,kalp,kact, 

*  xn,dxn,raa,orl,niiBr,BappartO,Bappart, 

*  a , t aat , a Igna , ataat , arror , arroznln , dalavnln , kllag , anda) 

ll(aba(arror).lt.apao)go  to  604 

11 (taat . la . 0 . 0) than 
llov«l 
datlsdalav 


l.davlstagt 

•llik>(f«Ta2-t«at)/(d«ni_2-d«l«v) 

g«VB_3«d«l«T 

ggnS^tast 

if (abs(slpk) .It. 0.02) than 

•lpk«(l«v*l-ta«t)/(dav«_l-dal«v) 

•ad  if 

d«l«Tsd«l«Y-t«st/slpk 

itcrontsitaroat't^l 

if (it«r2. cq. 1 .aad. itcroat .gt.30)th«n 
aostop-l 
go  to  692 

•ad  if 

call  dalav.atraaa (da .aoaop , apt . idia , aa . ipr apk , ir , i>t , daaO . 

-*■  itaroat .  iaaar ,  itar2 .  it  •r3 ,  idistl ,  zaro ,  val ,  cka ,  ckt ,  cs , 

+  iaaaz , apsi . apao. apaa, iral . aigaaO , iaiga.f aataraga ,f aav.car . 

dalav ,  at ,  aatra .  atra ,  dalzO  .dalyO ,  dalzO ,  dalz ,  dal  j ,  dalz ,  axac ,  aab , 
+  r.z.r.z.za.ya.za.zt.jt.zt.az.ay.az.aza.apa.aza.axt.ayt.azt. 

+  az,ay,az.aza,aya,vza.Bzt,ayt,8zt.aolidT,volaaa,daaaity, 
faO.fal.ftO.ftl.dft.dl.afz.afy.afz.aaz.aay.aaz.itrial, 

*  a,b.ak,bk,adjc,kxz,kyy,kzz,BX,ay,az.aap0.aap,aboB,kalp,kact, 
xa,dxa,raa,orf  ,aaar,aappartO,ai^part, 

*  a , t aat , a ig^a ,  ataat . arror , arrozaia , dalavaia , kf lag , aada) 

if (aba(arror).lt.apao)go  to  604 
go  to  929 

•laa 

ihigh^l 

davb^dalaT 

f.davhstaat 

alpk*(fata2-taat)/(davB_2-dalav) 

gava_3=dalav 

gavB3«taat 

dalavsdalav-taat/alpk 

itarout-itaroat-i-1 

if (itar2. aq. 1 .aad. itaroat .gt. 30)thaa 
aoatopsl 
go  to  692 

aad  if 

call  dalav.atraaa (da, aoaop, apt , idia,Ba, iprapk , ir , iat ,daaO , 

-*■  itaroat ,  iaaar ,  itar2 ,  itar3 ,  idiatl  ,zaro ,  ral ,  cka ,  ckt ,  ca , 
iaaaz ,  apa  i ,  apao ,  apaa ,  ir  al ,  a  igaaO ,  ia  iga ,  f  aavaraga ,  f  aav.car , 
dalav ,  av ,  aatra ,  atra  ,dalzO  ,dalyO  ,dalzO  ,dalz  ,daly  ,dalz  ,Bzac  ,anb , 
+  r,z,y,z,za,ya,ai,zt,yt,zt,nz,ay,iiz,aza,nya,aai,uzt,uyt,azt, 
BZ,vy,BZ,«za,vya,vza,vzt,«yt,vzt,aolidT,Toloaa,daaaity, 

+  f aO ,f al , f tO,f tl ,df t ,dl , af z, af y , af z, aaz , aay , aaz , itrial , 

a,b,ak,bk,adjc,kzz,kyy,kzz,az,ay,az,aapO,aap,aboa,kalp,kact, 
za,dza,raa,orf  ,aaar,aappartO,aappart, 


93 


*  •,t«at,>ipui,stMt,«rror,crrozBia,d«l«niia,kflaK,aada) 

if (abt(«rror).lt.«pao}go  to  604 
go  to  929 

and  if 

929  if  (taat.lo.O.O)tliaa 
iloasl 
davl-daloT 
f.doTlstaat 

if  (ilo«*iliigli.aq.l)go  to  940 

alidc=(gavB3-teat)/(gaTB_3-delav) 

gavB_3sdolav 

gam3staat 

if (aba(alpk) .lt.0.02)tban 

alpk- (f  av>2-taat )/ (davB_2-dalaT ) 
and  if 

dalav^dalav-taat/alpk 

itaroutsitaroQt-i-l 

if ( it ar2 . aq . 1 . and . it aront . gt . 30) than 
noatopsl 
go  to  692 

and  if 

call  dalav.atraaa  (dn  .nonop  ,npt ,  idin,  nn,  iprapk ,  ir ,  iat ,  danO , 
itaront ,  innar ,  itar2 ,  itar3 ,  idiatl  ,zaro ,  aal,  ckn,  ckt ,  ca , 

*  inmax , apa i , apao , apaa, iral , aigmaO , iaign , f navaraga , f nav.cnr , 

-«■  dalav .  av ,  aatrn ,  atrn  .dalxO.daljrO,  dalxO ,  dalx ,  daly ,  dalx ,  mxnc ,  nnb, 
r,x,y,z,xa,yn.zn,xt,yt,xt,ax,ny.nz,nxB,uym,nzm,nxt,iiyt,nxt, 

«x ,  ay  ,az ,  vxm,  vyn,«zm,«xt  ,«yt ,  azt ,  aolidv,  ▼olnma ,  danaity , 

fnO,fnl,ftO,ftl,dft.dl,afx,afy,afz,aaXtaay,anz,itrial, 

a,b,ak,bk,adjc,kxx,kyy,kzz,nx,ny,nz,BapO,aap,nbon,k8lp,kact, 

*  xn,dxn,ra8,orf ,nTimr,nappartO,Bappart. 

a ,  t  eat ,  a  igma ,  at  aat ,  error ,  arrormin ,  dalavmin ,  kf  lag ,  anda) 
if (ab8(error).lt.apao)go  to  604 
go  to  929 

elaa 

ibigksl 

davh=dalaT 

f.davhstaat 

if (ilo«aihigk.aq.l)go  to  940 
8lpk= (gevB3-taat)/ (gavB_3-dalax) 
gavn_3sdalav 
gaTB3=ta8t 

if (aba (alpk) .It . 0. 02}tban 

alpks(favm2-ta8t)/(davB_2-dalaT) 
and  if 

dalav^dalaa-taat/alpk 

itarontsitaront'fl 


if (itar2. aq. 1 .and. itaront .gt .30)than 
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noatop«l 
go  to  692 

•ad  if 

call  d«l«v_str«aa(dB,Maop,npt,idiB,aB, iprapk, ir , iat .danO, 

*  it aroat , iaaar , it •r2 , it  arS , idiat 1 , xaro . val , cka . ckt , ca , 

*  iBBaz , apai , apao.apaa, iral.aigaaO , isiga.fnaTaraga ,fnaT_cnr , 

*  dalaT,aY,Batra,atra.dalxO,daly0.dalzO,dalz,daly,dalz,axnc,nnb, 
+  r.z.y.z.XB.ja.za.xt.jt.zt.uz.oy.az.iuai.api.uza.uxt.ujrt.uzt. 

vz.vy.n.vxa.ayB.aza.axt.vyt.szt.solidT.Toliuia.danaity, 

*  fnO,fnl,ftO,ftl,dft,dl,afx,afy,afz,aax,aay,Baz.itrial, 

*  a,b,ak,bk,adjc.kxx,kyy,kzz,ax,ny,az,BapO,Bap,nbon,kalp,kact, 

*  xa,dxn,raa,orf .nnar.aappartO.Bappart. 

*  a ,  t ast ,  a ipaa ,  ataat ,  arror ,  arrozBin , dalavBin .  kf lag ,  aada) 

il(aba(arror).lt.apao)go  to  604 
go  to  929 

and  if 

and  if 

940  ddav=davh-davl 

•rroxBinslOO.O 
dalavBinsO.O 
itar2»2 

do  924  itarotttai,ontBax 

rtflap*davl+ddav*f_daYl/(f_daYl-f_daYh) 
if (aba (rtf lap-dalavt ) . It . 0 . 000000 1 ) than 
dalavsrtf lap*! . SO 

•laa 

dalav^rtflap 
and  if 

dalevt=dalav 

call  dalav_atra8a(dB,Bonop, npt.idiB.nB, iprapk, ir, iat, danO, 
itaront ,  inner ,  itar2 ,  itarS ,  idiatl ,  zero ,  ral ,  ckn ,  ckt ,  ca , 

*  inaiax , apai , apao , apaa, iral ,aigBaO, isigB,f navaraga ,fnav_cnr , 

*  dale V ,  ,  astm ,  atm , dalxO  ,dalyO ,  dalzO , dalx ,  daly ,  dalz ,  Bxnc ,  nnb , 

r,x,y,z,ZB,ya,ZB,xt,yt,zt,ux,ny,az,nxB,nyB,naB,uxt,uyt,uzt, 
vx,wy  ,vz,«xB,ayB,wzB,«xt,vyt,«zt,aolidY,volnBa,danaity, 
fnO,fnl,ftO,fti,dft,dl,afx,8fy,sfz,aBX,8By,8BZ,itrial, 

+  a,b,ak,bk,adjc,kxx,kyy,kzz,nz,ny,nz,BapO,Bap,nbon,k8lp,kact, 
xn,dxn,ra8,orf  ,nnBr,BappartO,Bappart, 

-i-  8 ,  t  •8t ,  8  i^a ,  ataat ,  error ,  arrozBin ,  dalavBin ,  kf  lag ,  aada) 
if (abs(arxor).lt.ap8o)go  to  604 

f_fsta8t 

if (f _f . It . 0 . 0} than 
davlsdalav 
f_davlaf_f 

•!•• 

davhsdalav 
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•ad  if 

dd«v3d«vh-daTl 
924  coatiau* 

itar3«it«r3-t'l 

it«roat«l 

it«r2»l 

if ( it 6x3 . It . 2} thaa 
go  to  692 

•la*  if ( it •r3 . aq . 2 . aad . isida . gt . 0) thaa 
iaagas'l 
go  to  692 

alaa  if (itar3 . •q.2 . aad.isida.lt . 0)thaa 
iaaatsl 
go  to  692 

also  if ( it •r3 . aq . 3 . aad . loaatra . It . 2 ) thaa 
go  to  692 

•Isa  if ( itar3 . aq . 3 . aad .highstra . It . 2) thaa 
go  to  692 

also 

vrita(99,*) ’sorry,  I  caaaot  fiad  a  satisfactory  dalav' 
vrita(99,*)'to  balaaca  tha  strass,  stop!  aad  of  outloop' 
stop 
aad  if 

c — updatiag  tha  positioa,  forca,  aappiag  ate.  for  tha  aaxt  stap 
c— conputiag  aad  storing 

604  call  vac_$eopy(xt,x,apt) 

call  vac_$copy(yt,y,npt) 
call  vac_|copy(zt,z,npt} 
call  vac_9icopy(Bap,aap0,npt} 
call  vac_$icopy(Bappart,BappartO,nxnc) 
call  matx_copy(apt,apt,fal,faO) 
call  rBatx3_copy (apt , apt , da , f t 1 , f tO ) 

dalxO=dalx 

dalyO=daly 

delzO=dalz 

dx»dalxOanB 

dysdalyOanB 

dz^alzOanB 

av=ar+dalav 

dalar_last=dalar 

if (dalav . gt . 0. 0) isida0=i 

if (abs (dalav) .It . 1 .0a-10)isida0s0 

if (dalav . It . 0 . 0} isidaOs-1 

if (iprapk .ga . 3)thaa 

if (abs (stra.c) . gt , 1 . 0a-8}slop=delav/stra_c 
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•ad  il 


c 

c 

c 

c 

c 

c 

c 

c 


call  v«c_|iiBit(kphi,12.0) 
call  v«c.$iiBit(kth«ta.24,0) 

il((ipr«pk.gt . 1 .and. lit .eq. (idi>t2« 
(ift/idiat2)}).or.kflag.flq. Dtkan 
call  riBit2(Bpt,Bzc.ctpz,z*ro) 
call  riBit2(Bpt,axc,ctpy,zaro) 
call  riBit2(npt,aacc,ctpz,zaro) 
call  iBit2(Bpt,Bxc,tkka,izaro) 
call  init2(npt,Bzc,thkt,izaro) 
and  if 


c  do  607  i*l,npt-l 

c  do  608  jsi-fl.npt 

c  dl(j.i)*dl(i.j) 

c  608  coat Inna 

c  607  continaa 


call  copjhalf (npt.dl) 

call  rinit2(dB,dB,nij ,zaro) 

call  rinit2(npt,npt,co&d.zaro) 
call  vac.linitCbrz, apt .zero) 
call  Tac_linit(br7,npt,zaro) 
call  vac_$iBit(brz,npt,zaro) 

call  init2(npt,npt,ncont,0) 
call  rinit2(npt,npt,tfor,zaro) 
call  rinit2(npt,npt,tft,zaro) 


condBaz=0.0 
condnin^ 100.0 
condiiaz2=0 . 0 
condnin2=100.0 
avr_fn=0.0 
navr_ln=0 

avr_cond2=0 . 0 
navr_cond2=0 


do  221  np=l,npt 
ih(np)=0 

do  222  j_kn=:l.tinb 

knanappart (adj  c (napCnp) , j_kn) ) 
if (kn.aq.O)go  to  222 

if (fnl(np,kn) .ga.-ap8a)go  to  222 


nij  ( 1 ,  l)anij  (1 ,  l)-fnx(np,kn)*Bx(np,kn) 
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Bij  (1 , 3)Bni  j  (1 ,2)-fBx(i^,lm)*By(Bp.kB) 
nij  (1 , 3)»nij(l  ,3)-t'iiz(np.kB)*nz(Bp.kn) 
ni  j  (  2 , 1 )  sBi  j  (2 , 1 ) -f  By  (up ,  kB)*nz  (np ,  ka) 
ni  j  (2 , 2)»Bi  j  (2 , 2)4-By  (19  ,kB)*By  (ap.ka) 

Bi  j  (2 . 3)»Bi  j  (2 , 3)-t-By  (19  .kB)*BZ  (bp  ,kB) 

Bij(3,  l)>Bij(3.  l)-»'Bz(i^,kB)*Bx(Bp.ka) 

Bij  (3 , 2)»Bi  j  (3 , 2)-(’Bz(bp  .kB}*By  (bp  ,kB} 
Bij(3,3)BBij(3,3)-t'Bz(i9,kB)*Bz(Bp,kB) 

■COBt(Bp,kB)si 

tfor(Bp,kB)ssqrt(lBl(Bp.kB)**2-t-ltl(ap,kB,l)**2-i- 
ltl(Bp,ka.2)*'*2-t-itl(Bp.kB.3)**2) 
il (abs ( ca ) .  gt . 1 . Oe-5 ) thea 

tit(ap,kB)saqrt(ltl(Bp,ka,l)«*2-»' 
ft  1  (ap ,  ka ,  2  )  a*2-»^lt  1  (ap .  ka ,  3  )  **2) 
aad  if 


kx-kxx (aap (ap) , j  Jta) 
ky=kyy(Bap(Bp) , j_kB) 
kz=kzz(Bap(Bp) , j_ka) 

xl=kx*Ba*delx-t-ky*taa(88tza(2 , 1)  )*nB*dely-t- 
kz*taa(88tra(3 . 1) ) ^aaadelz 

y  l*ky  *B«*d8ly+kx*taa(  aatra  (1,2))  *aB*delx-«- 
kz*taa(88tra(3 . 2) )*aa*delz 

zlskz*Ba*d8lz-*-kx8taa(88tra(l ,  3)  )*Bn*delx-i- 
ky*tan(88trB(2,3))«BB*dely 
xka-xt(ka)-xl 
ykB*yt(ka)-yl 
zka=zt(ka)-zl 

call  d8gra8_plii(Bz(Bp,ka).kpbi) 
call  degree_theta(ax(Bp,ka) ,By(ap,ka) ,ktheta) 

if  (  (  iprepk .  gt .  1 .  aad .  i8t .  eq.  (  idi8t24' 

(i8t/idi8t2)))  .or.kflag.eq.Dthea 
ih(ap)»ih(Bp)-i-l 

ctpx (ap , ih (ap) )=(r (ka) *xt (ap) +r (ap) *xkB) /dl (ap , ka) 
ctpy (ap , ih (ap) )=(r (ka) ♦yt (ap) +r (ap) *ykB) /dl (ap , ka) 
ctpz (ap , ih(Bp) )= (r (ka) *zt (ap) +r (ap) *zkB) /dl (ap , kn) 

ead  if 

if ( ipr «pk . «q . 3 . or . ( iprspk . la . 2 . aad . kf lag . 8q . 1 ) ) thaa 

fraal=ab8  (fal  (Bp,kB)abarRabarka) 
coad (ap , ka) =coBbaaxp ( coBkalog(f raal ) ) 

avr_coBd2=avr_coad2-i-coBd(ap,kB) 

aavr_coad2saavr_coad2'«-l 


axr.f B*avr_f a+f al(ap .ka) 
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navr_f  ii=n»vr_f  n+ 1 

c  if (kx.n«.0.or.k7.n«.0}cond(np.lm)s0.0 

if  (con<l(np ,  kn)  .  gt .  condaax)  condux=cond  (np  ,kn) 
if (cond(np , ka) . It . condain . and . cond(np , kn) . gt . 1 . 0«-20 ) 
-*■  condsin-cond(np,kn) 

brx(np) sbrx (np) -cond (np , kn)* (xkn-x (np) ) *barR*grad 
br y (np) =bry (np) - cond(np , kn)* (ykn-y (np) ) *barR*gr ad 
brz (np) »brz (np) -cond (np ,kn)* (zkn-z (np) ) *barR*grad 

and  if 

222  continue 

221  continue 


if (iprepk.eq.S.or. (iprepk.le.2.and.kflag.eq. l))then 
do  224  i=l,npt 

ptB=0 . 0 

do  223  j=l,npt 
ptB^ptB-i-cond  (  i ,  j  ) 

223  continue 
cond(i.i)=-pta 

224  continue 

call  conduct (npt,nxnc,nnb,dx,dy,dz,s8tm,x,y,z,r, 

+  nx, ny,nz, cond, brx.bry, brz, cur, voltf,c_eff, 

cxn ,  cr  es ,  cdxn ,  nap ,  nappart ,  ad  j  c  ,kxx ,  kyy ,  kzz ,  grad ,  barR) 

nr  it  e  ( 30 ,  *)  ist ,  sstm.c ,  condaax ,  condain 
write(30,*) ’  ' ,c_eff (1,1) ,c_eff (1,2) ,c_off (1 ,3) 

write(30,*) ’  * ,c_eff (2,1) ,c_eff (2,2) ,c_off (2,3) 

writ e (30 , *) '  * , c_ef f (3,1), c_ef f (3 , 2) , c.ef f (3,3) 

writo(30,*)’  ’ 

c -  coaputing  conductivity  by  nean  field  theory 

avr_f n=avr_f n/real (navr_f n) 

avr_f real=abs ( avr_f n*barR*b2urkn) 
avr_cond=conb*  exp ( conk*log (avr.f r eal ) ) 

cnc-2 . 0*r2*r2*avr_cond/(dx*dy*dz*barR) 

avr_cond2=avr_cond2/real(navr_cond2) 
cnc2=2 . 0*r2*r2*avr_cond2/ (dx*dy*dz*barR) 
write(30,*) '  ’ ,avr_fn,avr_cond,avr_cond2,cnc2*nij(3,3) 

c  uriteOO,*)’  '  ,avr_fn, avr_f real, avr_ cond, cnc 

vrite(30,*) •  ' , cnc*nij (1 , 1) ,cnc*ni j (1 ,2) , cnc*nij (1 , 3) 

vrite(30 , •) •  ' , cnc*ni j (2,1) ,cnc*ni j (2,2), cnc*ni j (2,3) 
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n:it«(30 , •) *  • , cac«nl j (3,1) ,c«c*iii j (3,2), c«c*ni j (3,3) 

»rit«(30,*)'  • 

end  if 


do  323  isl.da 
do  324  j’^i.da 

ni j ( i , j ) =ni j ( i , j ) / (2 . O*nbon) 
continue 
continue 


f nuxsO . 0 

fmin=10.0 

fnnax=0.0 

fiuiiin=10.0 

ftmax=0.0 

ftmin^lO.O 

do  258  i=l,npt-l 
do  259  j=i-fl,npt 

il(fnl(i, j) .gt.-ep8a)go  to  259 
if (tf or (i , j ) . gt . fnax) then 
fBaxstfor(i, j) 
inax-i 
j»ax*j 
end  if 

if (tfor(i, j) .lt.fBia)then 
fmin®tfor(i,j) 
imin^i 
jmin*j 
end  if 

if (ab8(fn0(i,j)) .gt.fBBax)fnnaz=ab8(fnO(i, j)) 
if (ab8 (fnO ( i , j ) ) . It . f nnin) f nBin=ab8 (f nO ( i , j ) ) 
if  (tf  t  (  i ,  j  ) .  gt .  ltBax)f  tB>ax=tf  t  (i ,  j  ) 
if (tft(i, j) .lt.ftiBin)ftBin=tft(i, j) 
continue 
continue 

if (ir . eq.  1 . or . kf lag . eq. l)then 

write (70 ,*) ir , i8t , nbon, fnmax , fnnin , -avr_fn 
do  358  i=l,npt-l 
do  359  j=i+l,npt 

if (fnl(i, j).gt.-ep8a)go  to  359 
«rite(70,*)fnl(i, j) 
continue 
continue 


il (ist . tq . idistl* ( iat/idist 1 ) . or . Mlag . eq. l)then 


«rita(54,*)ir,ist 
do  855  k=1.12 

lraction=kphi(k)/(2.0*nbon'fl.0a-10) 
vrita(54,863)k«15.0, fraction 
863  foniat(2x,i6>fl0.3) 

855  continue 

write (56, e)ir,i8t 
do  856  k=1.24 

f raction=ktheta(k)/(2 .Oenbon+l . Oe-10) 
vrite(56,863)k*15.0, fraction 

856  continue 


c  vrite(51,*)ir,ist 

c  do  589  i=l,npt 

c  vrite(51,57}uxt(i),uyt(i),uzt(i) 

c  Brite(61,57)*xt(i),wyt(i),wzt(i) 

c  57  format(2x.3fl5.6) 

c  589  continue 


c 

c 

c 

c 

c 

c  266 


c 

c  267 
c  265 


if (f nin . ge . epaa . and . f min . ne . 10 . 0) then 
ifrite(92,*)iat,  *  in,  ft* 
vrite(92,e)*  Smalleat  force: ' ,imin,jmin, 'fmins 

vrite(92,*)*  Largeat  force: * ,imax,j max, 'fmax- 

write(92,*)’  Smalleat  normal  force:*,*  fnmin^ 

vrite(92,e)*  Largest  normal  force: *, *  funax^ 


*  ,fmin 

*  ,fmax 

'  ,fxmin 

*  ,fnfflaz 


do  265  i=l,npt 
write(92,e}i 
do  267  j=l,npt 

if (fn0(i, j) .lt.-epsa)tben 

¥rite(92,266)  j  ,fnOU,  j)  ,ft0(i,  j,  1)  ,ft0(i,  j  ,2)  ,ft0(i,  j  ,3) 
format(5x,i4,4(f 15.11,lx)) 
end  if 
continue 
continue 


end  if 
end  if 

vr  it  e  (55 , 854)  ir ,  iat ,  8stm_c ,  nbon ,  kslp ,  kact 
854  format(2x,2iS,fi2.6,3i5) 

write(68 ,  *)-sstm_c  ,-fnavorage 
write(69,e)-sstm_c,-fnav_cur 

npaTt=0 
do  189  k=l,npt 

if (map(k) .ge.l.and.map(k) .le.ncell)then 
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npart^npart-t-l 
•ad  if 

189  continn* 


write (53 , 851 ) ir , iat .apart , tatra_c , dalev . 
d«asity,slop,«v 

851  forBat(lz,3i5.5fl0.5) 

nrit«(53,852)8(l,l),a(l,2),t(l,3) 
Brit«(53,852)a(2.1).a(2.2),8(2.3) 
vrit«(53.852)s(3,l).8(3.2),s(3.3) 
write(53,*)’  ’ 


Brit«(53,852)aij(l,l),aij(l,2),aij(l,3) 
«rit«(53.852}aij(2,l).aij(2.2),aij(2,3) 
write(53.8S2)aij (3,l),aij(3,2).aij(3,3) 
852  foraat(5z,3(lz,«15.7)) 


c 

c  + 

c 

c 

c  + 

c  + 

c 
c 

c  339 

c  489 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c  + 

c 

c 


if ( ( iprepk .gt . 1 . aad. iat . eq. (idiat2a 
(iat/idiat2) )) .or .kflag. aq. l)th«a 

write (80 .«) ir , iat . delzeaai , delyean , delz*im , 
aatra(l,2) ,aatra(l,3) ,aatra(2.3) 
aatra(2,l) ,aatza(3,l) ,aatra(3.2) 
do  489  i=l,apt 

write(80,339)i,r(i) ,z(i) ,y(i) ,2(i) 
fonBat(i5,f6.2,3(lz,fl0.5)) 
coatiaue 

write(80,*}kflag 

f adif s (f anaz-f Bmia}/4. 0 
ftdif=(ftBaz-ftBia)/4.0 

call  iait2(apt,BZc,thka,izero) 
call  iait2(Bpt,mzc,thkt,izero) 

do  702  i=l,apt 
jq=0 

do  701  j=l,apt 

if (fa0(i, j) .gt.-epaa)go  to  701 
jq=jq+l 

if (aba(fadif) .le.epaalthea 
thka(i, jq)sl 

elae 

thka  (  i ,  j  q)  -  iat  (  (aba  (f  aO  (  i ,  j  )  ) -f  aaia) /f  adif  )  ■•■  1 

•ad  if 


if  (ab8(tft(i,  jD.gt.epaa.aad. 
aba (f tdif ) . le . epaa) then 
thkt(i,jq)=l 

elae  if (aba(tft(i, j)) .gt.epaa.aad. 
aba (f tdif ) . gt . epaa) thea 
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c 

c 

c 

c 

c 

c  701 

c  702 

c 

c 

c 

c 

c  87 

c  86 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c  88 
c  89 
c  97 
c 
c 

e 

c 

c 

c 

c 

c 

c 
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tUrt  ( 1 ,  j  q)  >lnt  (  (»bs  (tit  (  i ,  j  )) -f  tain) /ltdil  )■«■  1 

•la* 

thkt(i,jq)*0 
•nd  il 

contiana 

contiana 

vrita(81,*)iat 
do  86  isl.apt 

vrita(81,87)i,ih(i) 

loraat(2x,2i5) 

coatiaaa 

«rita(82,a}i8t 

arit a (83 , * ) ist , lamax , f aaia , Itaax , Itmia 
do  97  i=l,apt 

*rita(82,88)(ctpx(i, j), j=l,10) 
vrita(82.88)(ctpy(i.j),j=1.10) 
nrita(82,88)(ctpz(i. j). jsl.io) 
vrit a (83 , 89) (thka(i ,j).j=l,10) 
vrita(83,89)(tbkt(i,j). j=l,10) 
lonuLt(10(lx,16.3)) 
lozaat ( 10( lx , 12) ) 
coatiaua 

aad  if 

if (iprapk . gt . 2)thaa 

clo8a(71 , 8tattt8='dalata  * ) 
clo8a (72 , 8tatU8=  *  dalat a ' ) 

opaa (aait=71,fila=' pradat 1 ' , atatua = ' aakaova  * ) 
opaa(aait=72, f ila^'coordl' , 8tata8= 'uakaovn ’ ) 
iar*l 
aad  if 

do  98  i=l,3 
do  98  j=l,3 
p88tra(i , j ) =88tra( i , j ) 
star ia( i , j ) =88tra( i , j ) 
if (iprapk. la. 2)staria(i,j)=0.0 
coatiaaa 

istol=i8t 

eTV=ax 

delav_la~dalav_la8t 
if (iprapk . la . 2)thaa 
istol^O 
avvsO . 0 
dalaT_la=0 . 0 
aad  if 
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if  (iprapk .  aq.  1 .  and.  abaCaitnO ,  ll+sgtraC  1 ,2)  ) .  It .  1 .  Oa-8 .  and . 
dans  it  7 .  ga .  (dantp-^ndan*daldan) .  and .  dana  ity .  It . 

(dansp-i-  (ndan-t’  1 )  *daldan)  )  tban 
ndan^ndan-i-l 
iwx*l 
and  if 

if (kflag.aq.l.or. (iat.aq.iatold+istpl.or.iar.aq. 1 .or. 

(icnan .  ga .  i_ciuni .  and .  iprapk .  aq.  1 )  )than 
i»r=0 

nrita(71 , *) ir ,nlO ,a20 ,n30 .pi .p2.dy , aolidv . signa 
arita(71,a)iat,iatol.atarin(l,l),atarin(l,2) ,starin(l,3) , 
8tarin(2 , 1 ) , star in(2 , 2) . atarin(2 , 3) . 
atarin(3,l),atarin(3,2),8tarin(3,3),aTv,dalaa_la,iaida0 


do  192  1=1, npt 
ntansQ 

do  191  k=l,npt 

ntaB=acont  (  i ,  k)  -fntaa 
continna 

nrita(71,*)i.nt«ai 
do  193  j=l,npt 

if(fnO(i. j).lt.0.0)tban 

»rita(71.*)j,fn0(i,j),ft0(i,j,l),ft0(i,j,2),lt0(i,j,3) 
and  if 
continna 
continna 

vrita(72,*)ir,ist,dx,d7,dz,nB,ijaO 
do  196  i=l,npt 

writa(72,a)i,r(i) ,x(i) ,y(i) ,z(i) 
continna 
and  if 

if (kflag.aq.l)go  to  699 
if (icnatn . ga . i.cnnn .and . iprapk. aq. 0)than 
if ( itr ial . ga . itrlala)tban 
vrite(99,a} ’■axlnua  nnnbar  of  trials  to  pack’ 

«rita(99,a) 'to  dasirad  density  is  raachad,  stop!’ 
stop 

alsa 

itrialsitrialtl 
go  to  199 
and  if 
and  if 

if (icnnn.ga.i.cnBn. and. iprapk. aq.l) go  to  599 
continna 

vrita(68,*)’ft’ 

«rita(69,*)’k’ 
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e  il(il>r«pk.«q.2)th«n 

c  »rit«(*,*)’irs’,ir,*,  d*Bsitjs> .density 

c  end  if 

699  contir.  .e 

stop 

snd 


c — subroutine  to  do  packing,  and  force  balancing 

subroutine  delev_stress(d>,noaop.npt , idia.nn, iprepk , ir , ist .denO, 
iterout ,  inner ,  iter2 ,  iterS ,  idistl  .zero,  val ,  ckn.  ckt .  cs . 

*  in>ax.epsi.epso.epsa.irel.sigB*0.isign.fnaTerage.fnav_cur. 

*  dele V .  ev .  sstm .  stm .  delzO  .delyO .  delzO .  delz .  dely .  delz ,  nzne ,  nnb . 

+  r.x.y.z.zn.yn.zn.zt.yt.zt.uz.ny.uz.uzB.uyB.uzB.uzt.uyt.uzt. 

sx.vy.sz.vn.vyn.sza.vxt.syt.szt.solidT.volune. density, 

+  fnO.fnl.ftO.ftl.dft.dl.sfx.sfy.sfz.snx.ssy.snz.itrial. 

-*■  a.b.ak.bk.adjc.kxx.kyy.kzz.nx.ny.nz.napO.nap.nbon.kslp.kact. 
xn.dxn.res.orf  .nunr.nappartO.aappart. 
s, test. sigsn.stest. error. erromin.delevnin.kf lag. astda) 

inplicit  real  (a-h.o-z) 
implicit  integer  (i-n) 

integer  dm 

real  x(npt),y(npt).z(npt).r(npt) 

real  xm(npt) .ym(npt) .zm(npt) .xtCnpt) .yt(npt) ,zt(npt) 
real  ux(npt) .uy (npt) ,uz(npt) .uxmCnpt) .nym(npt) .uzmCnpt) 
real  ux(npt) ,sy(npt) ,vz(npt) .vxmCnpt) .uyB(npt) ,Bzm(npt) 

Z€al  uxt(npt) .uytCnpt) ,uzt(npt).vxt(Bpt),vyt(npt) ,vzt(npt) 
real  a(idim, idim) ,b(idiB) .ak(8.6) .bk(6,6) ,xn(idim) 

real  f nO (npt . npt) . f nl (npt ,npt ) 
real  ftO(npt .npt.dm) .ftKnpt.npt.dm) ,dft(dB) 
real  sfx(npt} .sfy(npt} .sfz(npt).sBx(npt) ,sBy(npt) ,8Bz(npt) 
real  BxO,By0.az0.dl(npt.npt).dxB(idiB) 
real  nx(npt.npt) ,ny(npt,npt).nz(npt,npt).nxx,nyy,nzz 
real  res(idiB).s(dn,dB),s8tm(3.3).stm(3.3) 
c  real  r  t.7(idiB).s(dn,dB).sstm(6),8tm(6) 

integer  kxx(Bxnc.nnb) ,kyy(Bxnc.nnb) ,kzz(mxnc.nnb) 
integer  adjc(nxnt ,1."^) .Bap(npt) .BapO(npt) .stest 
integer  Bappart(Bxnc)  ,Bappaj:tO(Bxnc} 


ncell=nn*eda 

cknt=ckn-ckt 


c 


Brite(*.*) 'steps' .ist. '  outs' .iterout. '  iter2s' ,iter2. 
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c  ♦*  lt«r3s* ,it«r3 

c  «rit*(6,*)'  ioaar,  bz.  bj,  bz,  anab.  ratb,  nb,  aba,  atraaa' 

c  if ( iat . aq . idiat 1 * ( iat/idiat 1} ) tbaa 

e  *rita(49,*) 'atap»',iat, *  ont^’.itaront, *  itar2s ' , itar2 . 

c  ♦  *  itar3«' ,itar3 

c  arita(49,*)*  iaaar.  bz,  by.  bz,  aaab,  ratb,  nb, 

c  *  aba,  atraaa' 

c  and  it 

c - Bova  particlaa  according  to  naan  diaplacanaat 

factorzsatmd  ,1) 
lactorysatm(2,2) 
lactorzBatm(3,3) 

iactoradalav/da 

do  714  ial,npt 

iuai(l)3(lactorz-»^factor)az(i)'t’atzn(2,l)ay(i)-fgtzn(3,l)*z(i) 
nyB(l)>(factory»'factor)*y(i)'<-atzn(l,2)az(i)-t-atm(3,2)*z(i) 
nzn(i)s(factorz-«-lactor)«z(i)'«'atzn(l,3)*z(i)-i-atm(2,3)ay(i) 
«ZBt(  i)  satm(2  •  3) -atm(3 , 2) 
wya(i)satm(3,l)-atm(l,3) 
wzat(  i)  satm  ( 1 , 2) -atm  (2 , 1 ) 
zzdlazCD-t’Hzad) 
jad)syd)-t-ttyBd) 

ZB  d  )  «z  (  i ) ‘fuzz  d  ) 

ttz(i)B0.0 

vyCD^O.O 

uzd)s0.0 

wz(i)80.0 

¥yd)«0.0 

«z(i)s0.0 

714  continna 

dalz^dalzO*  ( 1 . 0-t'f  actorz-t’lactor) 
dalysdalyO* ( 1 . 0+lactory+lactor) 
dalz>dalzO«  ( 1 .  (Hlactorz-Mactor  ) 
danaitysaolidv/Cncalladalzadalyadalz) 

it dprapk . la . 1 . and . daaaity.lt . 0 . 5)th  i 
inBazoeinaaz/lO 

alaa 

inBazosinanz 
and  if 

c - innar  loop  to  balanca  tba  f orca 

ivaysl 

innarsQ 
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333  laatraimivr't'l 

call  v«c_|«dd_T«etor(za»az.apt,zt) 
call  ▼ac_tadd_vactor(jB,«y.Bpt,yt) 
call  T«c_ladd_vactor(xa,az,Bpt,st) 

call  v«c_tadd_T«ctor(tuai,iiz.Bpt,iixt) 
call  vac_|add_Tactor(nyB.ay,Bpt.nyt) 
call  Tac_tadd_vactor(iiai,«z,i9t.axt) 

call  vac_$add_Tactor(wxa,*z.Bpt,*xt) 
call  vac_tadd_T«ctor(wyB,*y,npt,«yt) 
call  ▼•c_|add_vactor(«XB,az«Bpt.«ct) 

11 ( iaaar . aq .  1 ) than 

call  Tac_$icopy(B^>0,M^,Bpt) 

call  vac_$icopy(BappartO,Bappart.wmc) 

alia 

call  uppingCapt ,aB,BZBC , aatra.dalx ,daly ,dalz , 
+  zt,yt,zt,Bap,Bappart, louts) 

If  (louts .  ga .  Dthaa 

«rlta(0O,*)*partlclas  fly  out  of  tha  systaa!!* 
nrlta(99,a)ir,lst.ltarout,iaBar,*  louts** .louts 
«r It a (99 . *) 1st , dalzaaa.dalyaaa , dalzaaa, 

*  sstrB(l.l).sstra(1.2),sstzB(l,3), 

*  sstra(2,l),sstra(2,2).sstrB(2,3), 

*  sstrB(3,l),sstrB(3,2},sstrn(3,3) 
do  334  1*1, apt 

*rlta(99,a)i,r(l) .xt(l) ,yt(l) .zt(l) 

334  coatlaua 

sad  If 
aad  If 

abf z*0 . 0 
abf y«0 . 0 
abf z*0 . 0 
abaa=0.0 
abaqr^O.O 
abBZ*0 . 0 
kslp*0 
kact*0 

call  rlalt2(apt,apt,fal,zsro) 
call  rlalt3(apt,apt,dB,ftl,zaro) 
call  rlBlt2(Bpt,Bpt,dl,Tal) 
call  ▼ac_tlalt(sfz,xq>t,zaro) 
call  ▼ac_$lalt(sfy,Bpt,zaro) 
call  vac_$lBlt(sfz,apt,zaro) 
call  vac_|lalt(sBz,Bpt,zaro) 
call  vac_$lalt(sBy,Bpt,zaro) 
call  Tac_$lBlt(sBa,apt,zaro} 
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if  (iniivr .  !• .  iimxo}th«ii 

call  riait2(idiB,ldia,a,saro) 
call  Tac_ti]iit(b,idia,xaro) 
aad  If 

aboasO 

if (ivaj.aq. I}faav«raga*0.0 
faaT.car^O.O 

c — coapatiag  graad  itiffacaa  aad  nabalaacad  forca  vector 
c - avalaatiag  coatact  forces  betseea  particle 

do  21  apsl.apt 
kap*6e(ap-l) 
do  22  j_ka«l,aab 

ka>Mppart(adje(Bap(ap) .  j_ka)) 
if  (ka.eq.O)go  to  22 

if (dl(ap.ka).lt.lOO.O.or.dl(ka,ap).lt.lOO.O)go  to  22 

kz>kzx(aap(ap) . 
kyak jy (aap (ap) , j  _ka) 
kxakzz(aap(ap) .  j  Jm) 

zlakz*aa*delz'*'k7*taa(sstra(2 , 1)  leaaedelp-t- 

*  kz*taa(setra(3,l))*aa*delz 

ylakyeaa*dely+kz*taB(sstra( 1 , 2) ) *a■edalz■^ 

*  kz*taa(sstra(3,2))*aa*delz 

zlakzeBa*delz-t'kz*taa(  sstra  (1,3))  *B»*delz■^ 
ky*taa(sstra(2,3))*aiii^dely 
zkaazt(ka)-zl 
yka«yt(ka)-yl 
zkaazt(ka)''zl 

dl2a (zka-zt (ap) ) **2+{yka-yt (ap) ) **2+ 

+  (zka-zt(ap))*e2 

if (dl2.gt. (r(ap)+r(ka)*0,001)*(r(ap)+r(ka)+0.001))go  to  22 
dl(ap,ka)asqrt(dl2) 

ovlapadl (ap , ka) -r (ap) -r (ka) 

if (ivay . eq. 2)thea 

ckasezp(aBda*log(abs  (ovlap/faaverage)  )  ) 
cktaO . 8*cka 
ead  if 


faacka*ovlap 
if (fa.gt . -epsa)tbea 
go  to  22 

ead  if 

fal(ap,ka)afa 


if  (ivaj .  cq .  l}f  aavaraf  •sfBaT«ra|;«-ff]i 
fnav_ciir«lBav_eiir4'fa 

nx (ap , kn) B ( xkB-xt (np) )/dl (np , kn) 
ay (ap . ka) « ( jkn-fe  (19) )/dl (np.kn) 
ax (ap , ka) « (zka-zt (ap) ) /dl (ap , ka) 
axx*ax(ap,ka) 

»y7*ayC»P»ka) 

azz«az(ap,ka) 

raxBr(ap)*axx 

raysr(ap)*ay7 

raz«r(ap)*azz 

axapBaxt(ap) 
ayiq>«ayt(ap) 
tizapsuzt(ap) 
vxapBX (ap) *«xt (ap) 
wyapBT (ap) a«yt (ap) 
vzapBx (ap) *«zt (ap) 

uxkaBnxt  (ka)  -  (f  actorx-i-f  actor)  *xl-a  tra  (2,l)*yl-ttrB(3,l)*zl 

ayka®ayt(ka)-(lactory+factor)*yl-atrB(l,2)*xl-atra(3,2)*2l 

azkasazt  (ka)-(factorz-»'f  actor)  *zl-stza  (1,3)  *xl-strB  (2,3)  *7! 

vxkBBr (ka) *vxt (ka) 

vykasr (ka) * vyt (ka) 

vzkaBr (ka) *azt (ka) 

nxrBuxka>uxap 

ayr=uyka-uyBp 

azr^uzka-ozap 

vxrBwxap4^vxka 

*yr*wyap+wyka 

«zr«vzap+vzka 

aboasaboa-*'! 

il(abs(c8) .lt.l.0«-8)thea 
ftl(ap,ka,l)=0.0 
ftl(ap,ka,2)B0.0 
ftl(ap,ka,3)=0.0 
kslpskalp-i-l 
also 

df t ( 1 ) =ckta ( ( 1 . -axx**2) •Tixr-Bxx*ayy*uyr-axx*azz*u2r 
-B2z*wyr+ayy*wzr) 

df t (2) ®ckt* (-Byy*Bxx*axr+ ( 1 . -ayy**2) *ayr-Byy*B22*a2r 
■•■azz*¥xr-Bxx*«zr) 

dl t (3) Bckt* (-azz*axx*Bxr-B2Z*nyy*ttyr+ ( 1 . -bzz**2) aazr 
-Byy*¥xr+Bxx*wyr) 

dft(l)Bdlt(l)->’ftO(ap,kB,l)*(l.-axx**2)-ftO(Bp,ka,2)* 
BXXfayy-ltO (ap , ka , 3) *Bxx*azz 


dft(2)sdlt(2)-ftO(Bp.k&,l)*Byy*iixzmO(Bp,kn,2)*(l. 

»yi*®yy ) -<to(Bp , kB , 3) •Bjy •Bz* 
dlt (3 ) «(Ut (3 }-f tO(Bp . ka , 1 ) *azz*azz-f to (ap , ka , 2) * 
Bzz*B  jy-ff  tO(Bp .  ka .  3)  *  ( 1 .  <-nzz*Bzz) 

lt««qrt  («t  ( 1)  **2-t-dlt  (2)**2-»dlt  (3)  **2) 

tz»dft(l)/ft 

ty»dit(2)/lt 

tz«dlt(3)/ft 

iKabsdt)  .lt.«paB)ftsO.O 

if (f t-cs*aba (fa) . gt . 0 . 0) thea 
ftscs*aba(fB) 
f t 1 (ap .ka , 1 ) sf t*tx 
f t 1 (ap . ka . 2)sf t*ty 
ftl(Bp,kB,3)sft*tz 
kalp^kalp-fl 
alaa 

ftl(Bp,ka,l)sdft(l) 
ftl(ap,kB,2)»dft(2) 
ftl(Bp,ka,3)=dft(3) 
kact=kact't'l 
and  if 

and  if 

fxsfnanxx+f tl (np.kn, i) 
f y»f n*nyy+f t 1 (np ,kn. 2) 
f  zsf  nanzz-»'f  1 1  (np  ,kn ,  3) 
if (aba(ca) .lt.l.0a-8)tban 
BXOsO.O 
ny0=0 . 0 
az0=0 . 0 
alaa 

BxOsrayaf t 1 (np,kn,3)-raz*f tl (np ,kn , 2) 
ayOsrazaf tl(np ,kn, l)-raxaf tl (np ,kn ,3) 
BzO=raxaf t 1 (np , kn , 2)-ray af t 1 (np , kn , 1 ) 
and  if 

8fx(np)=8fx(np)‘«'fz 
8fy(np)ssfy(np)^fy 
8fz(np)aafz(np)'t’fz 
aaut  (np)  aaax  (np)  ■•-axO 
any  (np)  saay  (np)  -ivy  0 
8Bz (np) ssBz (np) aazO 

af x(kn) =af X (kn) -f x 

8fy(kn)a8fy(kn)-fy 

af z (kn) saf z (kn) -f z 

88UC  (kn)  S8BX  (kn)  -huOar  (kn)  /r  (np) 

aay  (kn)  =8By  (kn)  ■•■ayOar  (kn)  /r  (np) 
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■■z(kB)  3>BZ  (kn)  ■»'azO*r  (ka) /r  (up) 

aMz^ablx-f  aba  (Iz) 

4bf j«abfy4abs(tj) 
abfzsablz'^abadz) 
abaac=al»Z'<'aba  (azO) 
abaiysabBj4^abs  (ayO) 
abaz^abaz-^aba  (azO) 

if (Innar.aq.  (iaaazo4^1))go  to  22 
kkn-6*(kn-l) 

ckat-ckn-ckt 

ak  ( 1 , 1)  =cknt*nzz*nzz-*-ckt 

ak(l  ,2)=ckiit*nzz*nyy 

ak  ( 1 , 3 ) =ckntaiizz*nzz 

ak(1.4)=0.0 

ak(  1 , 6)  a-ckt*iizz 

ak(l,6)<:ckt«nyy 

ak(2, l)=cknt*nyyanzz 
ak (2 , 2 ) = ( cknt*nyy*nyy +ckt ) 
ak (2 , 3) =ckntaoyy mzz 
ak(2,4)-ckt*xizz 
ak(2.B)-0.0 
ak(2,6)*-ckt*nzz 

ak(3 , 1 ) scknt*nzz*azz 
ak(3 , 2) =ckiit*szz*iiyy 
ak  (  3 , 3)  -  (  cknt*iizz*nzz't’ckt ) 
ak  (3 , 4)  =-ckt*iiyy 
ak(3,S)=ckt*nzz 
ak(3,6)=0.0 

ak  (4 , 1 )  =ckt*  (-ray*azz*BZZ'«’raz*nyy*nzz) 
ak(4,2)=-ckt*(ray*nzz*iiyy+raz*(l , 0-ayy*nyy) ) 
ak(4 , 3) =ckt* (ray* ( 1 . O-azzaazz)* raz*ayy*azz) 
ak (4 , 4) =-ckt* (ray*ayy+raz*azz) 
ak (4 , 6) =ckt*ray *azz 
ak (4 , 6 ) =ckt*raz*azz 

ak  (5 , 1 )  *ckt*  (raz*  ( 1 . 0-azz*azz)-«’raz*azz*azz) 
ak (5 , 2) =ckt* (-raz*azz*ayy*raz*azz*ayy ) 
ak(6,3)3-ckt*(raz*azz*azz-<-raz*(l  .0-azz*azz)  ) 
ak(5 , 4) =ckt*raz*ayy 
ak (5 , 6 ) 3-ckt* (raz*azz+raz*azz) 
ak (S , 6) =ckt*raz*ayy 

ak (6 , 1 ) =-ckt* (raz*ayy*azz+ray* ( 1 . 0-azz*azz) ) 
ak(6 , 2) =ckt* (raz*(l . 0-ayy*ayy)+ray*BZZ*ayy) 
ak  (6 , 3)  sckt*  (-raz*ayy*azz'<-ray*azz*azz) 


Ill 


ak(6 . 4)  «ckt*ruc*azz 
ak  (9 , 5)  sckt*ray*nzz 
ak  (6 , 6  )  »-ckt  *  (raz*iixx+ray  *ayy  ) 

knplsknp+l 

kap2skiip4^2 

knp3sknp*3 

kap4sknp<-4 

kiipS*knp4^5 

knp6«kjip+6 

kknlskkn-»'l 

kka2-kkn-f2 

kkaS-kkn-fS 

kkn4skkii4'4 

kknSskka+S 

kkn63kkBi-*'6 

a(knpl  ,knpl)>>a(kapl.kspl)-ak(l ,  1) 
aCknpl , knp2) »a(knpl ,knp2) -ak( 1 , 2) 
a(knpl  ,kBp3)sa(kiipl  .kiip3)-ak(  1 . 3) 
a(knpl ,  knp4)  «a(knpl  .]aip4)  ■t'ak(  1 . 4)  *r  (np) 
a(knpl .  knpS)  aaCknpl  .kspS) -t-akC  1 ,  S)  *r  (np) 
a(knpl .knp9)=a(knpl.knp6)'fak(l ,6)*r(np) 

a (knp2 . knpl ) sa(kap2, knpl ) -ak (2 , 1) 
a(kiip2 ,knp2)aa(kBp2,knp2)-ak(2 ,2) 
a(lmp2 .kap3) sa(knp2,kap3)-ak(2 , 3) 
a  (knp2 ,  kap4)  »a  (kap2 .  kap4) 'f  ak  ( 2 , 4)  *r  (sp) 
a(knp2 .  knpS)  3a(knp2  .kapS)  ■•■ak(2 ,  S)  *r  (ap) 
a(kBp2  .kap6) -a(kap2  ,kBp9)  ■•■ak(2 , 6)*r  (ap) 

a(kap3,kapl)=a(kap3,kapl)-ak(3,l) 
a(kap3 , kap2 ) sa(kap3 , kap2 ) -ak ( 3 , 2) 
a (]mp3 , kap3 ) sa (kap3 , kap3 ) -ak (3 , 3) 
a(kap3  ,kBp4}=a(kap3,kap4)-t-ak(3 ,4)*T(ap) 
a(kap3,kap6)sa(knp3,kap5)-fak(3,&)*r(ap) 
a(kBp3  ,kap6)sa(kBp3,kap6)-«-ak(3 , 6)*r  (ap) 

a (kap4 , kapl } sa(kap4 , kapl ) -ak(4 , 1 ) 
a(kap4 , kap2) sa(kap4 ,kBp2) -ak(4 , 2) 
a (kap4 ,  ]mp3 ) sa(kap4 , kap3 ) -ak (4 , 3) 
a(kap4  ,kap4)  =a(kap4  ,kap4}  ■<-ak(4 , 4)  ♦r  (ap) 
a(kap4,kBp6}=a(kap4,kap5)-»’ak(4,5)*r(ap} 
a(kap4 ,  kap6  )  sa(kap4 ,  kap6  )  ■•■ak  (4 , 6)  *r  (ap) 

a(kap6 .kapl ) -a(kap5.kBpl}-ak(5 , 1 } 
a(kap6 ,kap2)sa(kapS,kap2)-ak(6 , 2) 
a (kap6 , kap3 ) sa(kap5 , kap3} -ak (6 , 3) 
a(kap5 ,kap4) =a(kap5 . kap4)^^ak(6 , 4) *r (ap) 
a(kap6 , kap5) sa(kap5 ,kap5)  ■fakCS , 6)*r (ap) 
a(kapS ,  kap6  ) -a(kap5 ,  kap6)  ■•■ak  (S ,  6)  *r  (ap) 
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a (knpe . knp 1 ) *a(knp6 , knpl ) -ak C 6 , 1 ) 
a (lmp6 , knp2 ) «a ( knp6 , knp2 ) -ak ( 6 , 2 ) 
aCknpe , kap3) «a(knp6.knp3) -ak(6 . 3) 
admpe .  knp4}  «a(knpO ,  knp4)  4ak(6 , 4)  *r  (np) 
a(knp«  .knpS)  sa(ksp6  ,knp5) -f  ak(6 . 5)  ar  (up) 
a(knp6 , knpe) sa(k&p6 ,knp6) aakCe . 6) ar (np) 

a(lcnpl ,  kknl )  aak(  1 , 1 ) 
a(knpl,kkn2)sak(l,2) 
a(knpl,kkii3)sak(l,3) 
a(knpl .kkn4)=ak(1.4)ar(kn) 
aCknpl .kknS) aak(l ,S)ar(kii) 
a (knpl , kkn6) =ak ( 1 , 6) ar (kn) 

a(knp2.kknl)aak(2.1) 
a(knp2,kkn2)>ak(2,2) 
a(knp2,kkii3)-ak(2.3) 
a(knp2,kkn4)=ak(2,4)ar(k&) 
a (knp2 , kkn5} -ak (2. S)ar(kn) 
a(knp2 . kkn6) =ak (2 , 6) ar (kn) 

a(knp3,kknl)sak(3, 1) 
a(knp3.kkn2)=ak(3.2) 
a  Omp3 ,  kka3  )  =ak  (  3 , 3) 
a(knp3 , kkn4) =ak (3 , 4) ar (kn) 
a(knp3 , kknS) sak (3 . S) ar (kn) 
a (knp3 , kkn6) sak (3 , 6) ar (kn) 

a(knp4,kknl)3ak(4,l) 
a(knp4.kkn2)aak(4,2) 
a (knp4 , kkn3 ) =ak (4 , 3) 
a(knp4 , kkn4) =ak (4, 4) ar (kn) 
a (knp4 . kknS) =ak (4 , S) ar (kn) 
a(knp4 , kkn6 ) sak(4 , 6) ar (kn) 

a(knp5,kknl)=ak(5,l) 
a(knp6,kkn2)sak(S.2) 
a(knpS , kkn3 ) sak ( 5 , 3) 
a ( knp5 , kkn4 ) =ak ( 5 . 4) ar (kn ) 
a(knp5,kkn6)=ak(S.6)ar(kn) 
a(knp5 , kkn6) =ak (5 . 6) ar (kn) 

a(knp6,kknl)=ak(6,l) 
a (knp6 , kkn2 ) =ak (6 , 2) 
a (knp6 , kkn3) =ak (6, 3) 
a (knp6 , kkn4) =ak (6 . 4) ar (kn) 
a (knp6 . kknS ) -ak (6 , S) ar (kn) 
a(knp6,kkn6)sak(6,6)ar(kn) 

bk(l,l)=ak(l,l) 

bk(l,2)=ak(1.2) 

bk(l,3)aak(l,3) 
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bk(1.4)«-«k(1.4) 

bk(l,5)*-*k(l.B) 

bk(1.6)»-ak(1.6) 

bk(2.1)»ak(2,l) 

bk(2,2)*ak(2.2) 

bk(2.3)sak(2,3) 

bk(2.4)*-«k(2,4) 

bk(2,6)»-ak(2.5) 

bk(2.6)«-ak(2.6) 

bk(3.1)=ak(3,l) 

bk(3,2)=ak(3.2) 

bk(3.3)*ak(3,3) 

bk(3,4)*-ak(3.4) 

bk(3,S)a-ak(3.6) 

bk(3,6)=-ak(3.6) 

bk(4.1)*-ak(4,l)*r(ka)/r(np) 
bk (4 . 2) =-ak (4 . 2 ) *r (ka)/r (np) 
bk (4 . 3) =-ak (4 . 3) *r (ka)/r (np) 
bk(4 , 4)  =ak(4 ,4)  *r(kn)/r(iip) 
bk(4.S)=ak(4.S)*r(kn)/r(np) 
bk(4 , 6) s&k(4 ,6) *r (kn)/r(np) 

bk (S , 1 ) s-ak ( E , 1 ) *r (ka) /r (up) 

bk(S.2)s-ak(6,2)*r(ka)/r(np) 

bk(5.3)»-ak(5.3)*r(ka)/r(np) 

bk(5,4)*ak(S.4)*r(kn)/r(np) 

bk(6,S)*ak(S,8)*rOai)/r(np) 

bk(8,6)3ak(5,6)*r(kn)/r(np) 

bk (6 , 1 ) s-ak ( 6 , 1 ) *r (ka) /r (np) 

bk(6.2)=-ak(6,2)*r(ka)/r(iip) 

bk(6.3)=-ak(6,3)*r(ka)/r(np) 

bk(6.4)»ak(6,4)*r(kn)/r(np) 

bk(6 , 8)*ak(6 .8)*r(ka)/r(np) 

bk(6,«)=ak(6,6)*r(ka)/r(np) 

aCkknl ,kknl)=a(kknl ,kkal)-bk(l , 1) 
a(kknl ,kkB2) =a(kkal ,kka2}-bk( 1 , 2) 
a(kknl , kkn3)=a(kknl ,kkn3) -bk( 1 , 3) 
a(kknl  ,kkn4)-a(kknl  ,kk&4)-^bk(  1 ,4)*r(kii) 
a(kkal ,k)ai8)=a(kkal,kkn8)+bk(l ,8)*r(kn) 
a(kknl ,kkn6)=a(kknl,kkn6)+bk(l ,6)*r(kn) 

a(kkn2 , kknl ) 3a(kkn2 , kknl ) -bk(2 , 1 ) 
a(kkn2 ,kkn2)=a(kkn2,kkn2)-bk(2 ,2) 
a(kkn2 ,kkn3) sa(kka2 ,kkn3)-bk(2 , 3} 
a  (kkn2 ,  kka4}  =a  (kkn2 ,  kkn4) 't’bk  (  2 , 4)  ar  (kn) 
a(kkn2  .kkn8)»a(kkn2.kkn8)-i-bk(2 ,  E)*r  (kn) 
a(kkn2  ,kkn6)=a(kkn2,kkn6)-«'bk(2 , 6)ar  (kn) 
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a(k]m3  .kknl)3a(klai3.klml)-bk(3 , 1) 
a(kkn3 .kkn2) »a(klm3 ,kkn2)-bk(3 . 2) 
a(kkn3  ,kk]i3)  sa(klai3.kkn3) -bk(3 , 3) 
a(kkn3  .klm4}«a(kkB3.k]m4)-*^bk(3 .4}*r(lm) 
a(kkii3  .kknS)  «a(kla3  ,kkB5)  4^bk(3 , 6)*r  (kn) 
a(lckn3.kkn0)sa(kk&3.k)m6)-t^bk(3.6)*r(kii} 

a(klm4 , kknl } sa(klm4. kknl) -bk(4 . 1 ) 
a(kkn4 .  kkn2)  =a(kkii4 .  kk]i2) -bk(4 , 2) 
a  (]ckn4 ,  kkii3  )  =a(klm4 .  kkn3) -bk  (4 , 3) 
a(kkn4 .  kkii4  )  sa(kka4 .  kkn4) -^bk  (4 . 4)  *r  (kn) 
a(kka4  .kknS)Ba(kkB4.kknB)  4^bk(4 , 6)  *r  (kn) 
a(kkii4  .kka6)  «a(kkii4,kkn6)4^bk(4 . 6)*r  (kn) 

a(kkaS,kknl)>:a(kkiiE.kknl)-bk(5.1) 
a(kkn5 .kkii2)sa(kkB6.kkn2)-bk(5 .2) 
a(kkn6 , kkB3) sa(kkn6 , kkn3) -bk(6 , 3) 
a(kkii5 .  klm4)  =a(kkn5  .kkn4) -i-bkCS ,  4)  *r  (kn) 
a(kknS,kkn5)-a(kknS.kkii5)-t^bk(6.S)*r(kn) 
a(kknS ,  kkn6)  sa(kk&S  ,kkn6) '<-bk(6 , 6)  *r  (ka) 

a(kkn6  .kkiil)sa(kkn6.kkiil)-bk(6 . 1) 
a(klm6 .  klm2)  sa(kkn6.  kk2i2)-bk(6 , 2) 
a(kkn6 , kkn3) sa(kkn6 .klm3) -bk(6 , 3) 
a(kkn6 , kkn4) >a(kkn6 .kkn4)  -MdcCC . 4) *x (kn) 
a(kkn6 , kknS) »a(klm6,kkn5) *bk(6 . 6) *r (kn) 
a(kkn6  ,kkii6)sa(kkn6,kkn8}-»-Ut(6 , 6}*r  (kn) 

a(kknl ,knpl)=bk(l , 1) 
a(kknl.knp2)=bk(1.2) 
a (kkn 1 , knp3 ) sbk ( 1 , 3) 
a(kknl ,knp4)=bk(i,4)*r(np) 
a(kknl ,knp5)sbk(i,S)*r(np) 
a(kknl ,knp6)=bk(i,6)*r(np) 

a(kkn2 ,knpl)sbk(2, 1) 
a(kkn2 , knp2) =bk (2 , 2) 
a(kkn2 , knp3) =bk(2 , 3} 
a(kkn2 , knp4) =bk (2 , 4) *r(np) 
a(kkn2 , knpS) =bk(2 , S) *r(np) 
a(kkn2 , knp6} =bk(2 , 6) *r (np) 

a(kkn3.knpl)=bk(3.1) 
a (kkn3 , knp2} =bk (3 , 2) 
a (kkn3 , knp3) =bk (3 , 3) 
a(kkn3 , knp4) =bk (3 , 4} *r (np) 
a(kkn3 , knp5) =bk (3 , 5) *r(np) 
a(kkn3 , knp6) =bk (3 , 6) (np) 

a(kkn4,knpl)=bk(4, 1) 
a(kkn4 ,knp2) =bk(4, 2) 
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«(kkn4 , knp3) «bk (4, 3) 

&(kkn4 ,  lmp4) sbk (4 , 4) *r (np) 
a (kkn4 . knpS ) sbk(4 , 6) *r (np) 
a  (kkii4 ,  knp6)  sbk  (4 . 6)  *r  (np) 

a(klm6,knpl}3bk(S,  1) 
a(kknS.knp3)-bk(S.2) 
a(klm5.knp3)>bk(6,3) 
a(kkn5 . kBp4) -bk (6 ,4) *r (ap) 
a (kknS , knpS) sbk ( 6 , S) *r (np) 
a (kkaS , knp6 ) sbk ( S , 6) *r (np) 

a(kka6 , knpl ) sbk (6 , 1) 
a(kkn6 . knpS) -bk (6 . 2) 
a (kka6 . knp3 ) sbk (6 , 3) 
a(kkn6 ,knp4) -bk(6 ,4) *r(ap) 
a(kkn6 , knp5) -bk (6 . S) ar (ap) 
a(kka6 , kap6) =bk(6 , 6)ar(ap) 

b (kapl ) sb (kapl ) -f z 
b (kap2) sb (kap2) -ly 
b(kBp3) sb(kap3) -Iz 
b (kap4) -b (kap4) -azO 
b (kapS) -b(kap5) -ayO 
b(kap6) =b(kap6) -BzO 

b(kkal)»b(kkal)+lz 
b(kka2)»b(kka2)+fy 
b(kka3)sb(kka3)alz 
b (kka4) «b (kka4) -nzOar (ka) /r (ap) 
b (kkaS) -b (kka5)-By0ar (ka) /r (ap) 
b(kka6) =b(kka6) -BzOar (ka)/r (ap) 

22  coatiaue 

21  coatiaua 

sii]Bb=0.0 
do  804  i=l,apt 
iitp=6*(i-l) 

saab=8UBbaab8(b(iitp-»-l)  )-i-ab8(b(iitp-f2)  )+ab8  (b(iitp-f3)  ) 
804  coatiaue 

il(iaaer.eq. l)thea 
8unbl-8UBib 
ratbsl.O 

p8UBb5=8UBb 

ead  if 

if (mod ( iaaer , 5) . eq . 0)  thea 

ratbsab8 (p8UBb5-8UBb) / ouabl 
p8UBb5=8unb 
ead  if 
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sslzsvoc_#as«a(sf X ,npt) 
ssf  y=T«c_$uuB(af  7  .apt) 
aaf  zsv«e.tuaB(sf  z .  apt) 
aaaxsvac.laaaaCaaz .apt) 
aaBy«vac_taaia(aB7 .apt) 
aaBZBvac.laaaaCaaiz ,  apt) 

rtx=0 . 0 
rty=0 . 0 
rtzsO . 0 
rtBuc=0 . 0 
rtay^O . 0 
rtaz-O . 0 

if (ablz . aa . 0 . 0)rtz«aaf z/ (2 . 0*abf z) 
if (abf y . aa . 0 . 0)rtysaaf y/(2 . Oaabf y ) 
if (abf z .aa . 0 . 0)rtz=aaf z/(2 .0*abf z) 
if  (abauc  .aa .  0 . 0)rtuaaaaac/(2.  Oaabatz) 
if  (abay .  aa .  0 . 0)z-tBysaaBy/ (2 .  Oaabaiy  ) 
if (abaa . aa . 0 . 0)rtBUEsaaBz/(2 . 0*abaa) 

voltuaasncalladalzadalyadalz 
if (aod(iaaar,6) .aq.O)thaa 

call  atraaaCda.apt.aoaop.apaa.ax.ay.az.fal.ftl.r 
,  voliuia ,  a ,  a  igaa .  0 .  iaigw) 
aad  if 

*r it a (6 . 48 ) iaaar , rtz . rty , rtz , aaab . ratb . aboB , kalp > a igaa 
if ( iat . aq . idiat 1* ( iat/idiat 1) ) thaa 

vrit a (49 , 48) iaaar . rtz . rty . rtz , aaab . ratb , aboa .kalp . a igaa 
forBat(5z,i5.3f8.3.2fl0.5.2i5,f 12.8) 
aad  if 

lfl=0 

if ( ipr apk . la . 1 ) thaa 

if (iaaar .ga . 10 . aad . aigaa.lt . aigaaO) If 1=1 
if  (  iaaar .  ga .  10 .  aad .  (rtz+rty-t^rtz) .  It .  0 . 03)  If  1=  1 
alaa 

if (ca.Ba.O.O)thaB 

if (aaab . It . apa i*0 . 5 . or . ratb. It . apaiaO . 6 . or . (rtx+rty+rtz+ 
rtaz-t-rtayartaz)  .lt.0.42)lfl=l 

alaa 

if  (aaab .  It .  apai .  or .  ratb .  It .  apai.  or .  (rtz+rty-t-rtz)  .  It . 
0.21)lfl=l 
aad  if 

if (aigaa. It . 0 . SaaigaaO . and. iaaar . ga . 10)lf 1=1 
aad  if 

if ( iaaar . gt . iaaazo) If 1=1 

if (Ifl.aq.O)thaa 
go  to  666 
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it (ivay . aq. 2 . or . «L ■ (uMla) . It . 0 . 00001 )th«a 
ekii«l .  0 
cktaO.S 

t&aT_cTir«taaT_ciir/roal(&boB) 

c  «rito(*,*)'laavara£«*,taaT«ra(o 

c  writaCo.ol'taar.ciir’.laav.cor 

go  to  888 

•laa 

ivaj«2 

taavaragastaavaraga/raal  (aboa) 
taav_car*taav.cTiT/raal(aboa) 

iaaax'l 
go  to  333 
aad  it 

aad  it 

666  it(iprapk.gt.l.aad.sigBa.lt.0.6asigBa0.aad.aod(iaaor,S).aq.0)t]iaa 
go  to  814 

aad  it 

it (iprapk . la . 1 . aad.daasity.lt . 0 . S)tliaa 
apsdosapsialO.O 
iralosiralaO.S 

alsa 

apsdo«apsi 

iralosiral 

aad  it 

c— solviag  tha  qaasi-liaaar  systaa  ot  aqaatioas 

call  ralaxddia.apt ,iralo,apsdo,ort,a,b,za,aaBr,ras, 
dxa ,  ist ,  itarout ,  iaaar) 

c - updatiag  tha  tluctuatioa  displacaaaat 

do  382  itasl.apt 

ktB=3a  (da- 1 )  *  (  itst- 1 ) 
uz  (  ita)  =iiz  ( ita)  ■•'za  (kta+l ) 
uy  (ita)  suy  (  ita)  ■•■za(kta4>2) 
uz  (ita)  snz  (  ita)-i-za(kta4-3) 
vz  (  ita)  Bvz  (  ita)  ■*'za(ktaf4) 

«y  (ita)  =vy  (ita) -i-za  (ktaf-S) 
vz  (  ita)  =vz  (  ita)  ■l■za(ktBf■6) 

382  coavtnfjo 

go  to  333 

c - oaca  torca  balaaca  is  achiavad,  coaputiag  tha  strass 
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call  ■tr«si(dB,apt,Bonop.«psa.nz,B]r.BZ,fnl,ftl,r 
.Tolwa.a.iigMtltiaigB) 

c  il ( iat . aq . idiat 1* ( iat/idiat 1 ) )thaa 

c  «rita(49,47)iaaar,rtz,rty,rtz,rtBx.rtBj,rtBZ,aboa,kalp,algaA 

c  47  fonutdz, 14, 618. 2. 315.110.7} 

c  aad  11 

c  «rlta(08,649)lat,ltarout,laaar.aaBr,rtz,rt7,rtz,rtaz,rtay,rtaz 

c  649  loz«Lt(3z,31S,17.618.3) 

814  taat>algBia-algBaO 

11 (t aat . It . 0 . 0)thaa 
ataat»-l 
alaa 

ataatal 
aad  11 

arrorst  aat/a IgaaO 

11 ( It  ar2 . aq . 2 . aad . aba ( arror ) .lt.arTozHla)thaa 
dalavala^dalav 
arronlaaaba  (arror  ) 
aad  11 

c  vrlta(50,601)lr,lat,ltaroat,ataat,atra(3,3) , 
e  *  arror, dalav.av.daaalty 

11  (Iprapk .  la .  Dthaa 

wrlta<*,*)'ira',ir,',  trial* *,ltrlal,»,  lat«',lat, 

*  *,  daaalty*' .daaalty 

alaa  11 ( Iprapk. aq.2)thaa 

c  vrlta(a,a)*l  aa  ralazlag  tha  atraaa,  plaaaa  salt.* 

aad  11 

ar It a (a , 601 ) Ir , lat , Itaront , ataat , taat , arror , dalav , av , daaalty , la IgB 
601  lonBat(13.15.213,112.9.110.6.113.9.2110.6,15) 

c - chock  11  packlag  criteria  ara  aat  or  aot 

11 ( Iprapk . la . 1 . aad . daaalty . gt . daaO . aad . arror . la . apao . 

-i-  aBd.aba(aatra(3,l)-i-aatra(l,2)).lt.l.0a-8)kllag*l 

return 

end 

c - aapplag  partlclaa  lato  alcrocall  coordinate  ayatoa 

aubrontlaa  aapplng(apt,na,axnc,aatm,dalz,daly,dalz,x,y,z,aap, 
aappart ,  lonta) 
lapllclt  real  (a-h,o-z) 
lapllclt  Integer  (1-n) 

real  aatza(3 . 3) , z(npt) , y (apt) , z(npt) , dalz ,daly , dalz , dz ,dy , dz 
Integer  aap(npt),aappart(axnc),nB 


ds«r (bb) •d«lz 
dj«rMl  (aa)  *d«l  j 
ds«r sal (bb) *d«lz 
ioutsaO 

call  Tac.tiiBitCBi^part.azBCiO) 
do  1  isl.Bpt 

ZBia«taB(BBtzB(2 , 1) )*j(i)'ftaB(astra(3, 1) )*z(i) 
if  (z  (  i  } .  go .  mio)  thaa 

BZ*iat ( (z ( i ) -ZBia) /dz) 

also 

az>iat ( (z ( 1) -ZBia) /dz ) •! 
aad  if 

jBia>taa(sstza(1.2))*z(i)-i'taa(astra(3,2))ax(i) 
if  (y(i)  .8a.yBia)tliaa 

By«iat ( (y (i)-yBia)/dy) 

alaa 

By»iat((y(i)-yBia)/dy)-l 
aad  if 

zaia«taa( astZB ( 1 , 3) ) *z ( i ) > taaCaatra ( 2 , 3) ) *y ( i ) 
if (z(i) .ga.zaia)thaa 
az«iBt ( (z ( i ) >zaia) /dz ) 

alaa 

az«iat ( (z ( i ) -zain) /dz ) - 1 
and  if 

z(i)Bz(i)  -Bz*dz-By*tan(  astni(2 , 1}  }  *dy-BZ*taB  (  astrn  (3 , 1 )  )  adz 
y  (  i)  ay  (  i) -Byady-BZ*tan(aatrn(  1 , 2)  )  •dz-BZ*tan(aatm(3 , 2)  )  *dz 
z  (  i)  «z  (  i) -BZ*dz-BzataB(s8tzn(  i ,  3)  }  *dz-By*tan(  B8trn(2 , 3)  )  *dy 

zain3tan(s8tzn(2,l))*y(i)-«'tan(satrB(3,l))*z(i) 
yainataaC a atzn ( 1 , 2) ) *z( i ) atanC astrn (3 , 2) ) *z ( i ) 
ZBinataB(aatzn(l,3))az(i)-»’taB(a8trB(2,3))*y(i) 

kbaiat ( (z ( i) -zain) /dalz) a 1 
if  (kb .  aq.  BB-f  l)kb>nB 
jb*iBt(  (y(i)-yBin)/daly)-t’i 
if  (jb.aq.aBal)jbsaB 
ib^iat ( (z ( i) -ZBia) /dalz) ♦ 1 
if  (ib .  aq.na-f  1)  ib^nB 

ibar  >  iba  (  j  b- 1)  aaaa  (kb-1 )  *nB*nB 

if (aba (bz) . ga . 2 . and . aba (ay) .ga . 2. and . aba (bz) . ga . 2) than 
«rita(*,a)*particla  ',i,'  ia  locatad  ontaida* 
«rita(a,a)*kb,jb,ib,nB,ibar«  '.kb.jb.ib.na.ibar 
iontasionta^l 
and  if 


1 


■ap(i)«ibar 
■appart (ibar)«l 
eontiana 

ratara 

aad 


c - coaipatiag  aierocall  adjacaaey  aatriz  aad  ”kzz,kyy,kzx" 

aabrontiaa  adj3(Bxac,ijs,aa.nBb,adjc,kxz,kyy,kzz) 
iaplicit  raal  (a-b,o-z) 
iaplieit  iatagar  (1-a) 

iatagar  adjeCazac.aab) 

iatagar  kzzCazac.Bab) ,kyy(Bxac,nBb).kzz(Bzae,aab) 

do  60  k«l,aB 
do  60  j*l.aB 
do  40  isl.BB 

kl»i+ ( j -1 ) *Ba+ (k-1 ) *Ba*n* 
iadaxsO 

do  31  kk«k-lja,k-t'ija 
i<(kk.lt.l)thaa 
kz«l 

alaa  if (kk.gt.Ba)tbaa 
kz«-l 

alaa 

kz>0 
aad  if 

kkbarskk-*-kz*BB 

do  36  kj»j-ija,j+ija 
if(kj.lt.l)tbaa 
ky»l 

also  if(kj.gt.aB)thaa 
ky»-l 

alaa 

ky*0 
aad  if 

kjbar®kj+ky*a« 

do  30  kiai-ijs^i-fija 
ifCki.lt.Dthaa 
kzsl 

alaa  if (ki.gt.BB}thaa 
kx»-l 

alaa 

kzsO 
aad  if 

kibarski-t’kzana 
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k2«kib«r»  (k  jbur- 1 )  *aa-»  (kkbar- 1 ) 
if (kl.«q.k3)go  to  30 
iadox«iadoz-*'l 

a4jc(kl , iad«x)«k2 
kxx(kl , iBdox)«kx 
kyjCkl .  iitdox)«k7 
kzx(kl , iad«x)«kz 
coatlniio 
contiano 
coatiauo 
coatlaao 
coatiaao 
coatiaao 

rotaza 

•ad 


c - llasar  aqaatioa  aolvar 

subrontlaa  r«lax(idia.apr.lr«l,«p«d,orf .a,b,xa,ar,r«i,dxa, 
+  i*t , iont , iaa) 

implicit  real  (a-h,o-z) 
implicit  iatagar  (i-a) 

a(idia,idia)  ,b(idia)  ,raa(idim),dxB(idia)  ,za(idim) 
zaroBO.O 

call  ▼ac_liait(xatidia,zaro) 
call  vac_tiait(dxa,idia,zaro) 
call  vac_tialt(rai,idia,zaro) 

iaumsO 

100  iaamsiaam'*’! 

raax«0.0 
iaax«0 
rasamBO . 0 
if ( iaam . gt . 1 ) tbaa 
do  22  isi.idim 

if(abaU(i.jaax)).lt.l.0a-10)go  to  21 
ras  (i)  sras  (  D-^aC  i ,  jmax)  *dxa(  jmax) 

21  raaaa«raaaa+aba(raa(i)) 

if (abs(ras(i)) .gt.raax)tbaa 
raax«abs(ras(i) ) 
imax«i 
•ad  if 

22  coatiaaa 

•Isa 

do  20  iBl.idim 

if(abs(b(i)).lt.l.0a-12)go  to  20 


30 
35 

31 
40 
50 
60 
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r«i(i)>-b(i) 
r  •siiM-t'abs  (r«s  (  i)  ) 
if (abs(r«s(i)) .gt.raax)tbaa 
raax«abs(raa(i) ) 
iBax>i 
and  if 
20  continua 

raanBl«raaiui 

praaiui60«raaai 

and  if 

jaax«iaax 

dxii(  jaax)  «-raa  (  iaax) /a(  iaax,  jsax) 
z]i(jBaz)«zn(jBax)4dxn(juLa) 

if (aod ( inuB , SO ) . aq . 0 ) than 

ratiSOsaba (praanaSO-rastuil/rasial 
prasimSOaraaiai 

if (ratiSO.gt.apsd. and. iiiiia.lt. iraDgo  to  lOO 
Br*iiinB 
raturn 
and  if 
go  to  100 

aad 


c — shaffliag  algoritba 

subroutiaa  shaf f la(acall .ashuf .Iflip.arif , iaaad ,b,b1 ,kl ,k2} 

iaplicit  raal  (a-h.o-z) 

iaplicit  iatagar  (i-a) 

iatagar  B(acall) ,Bl(acall).kl(acall) ,k2(acall) 

do  65  isl,acall 
a(i)ai 
Bl(i)sO 
kl(i)=0 
k2(i)s0 

65  coatiaua 

acalsBcall/2 

aca2sBcall-acal 

do  100  jiBsl.aahaf 

do  105  If »1, If lip 

aids iat (raa( is  aad) aacall) 
do  102  isl.acall 
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if (i . !• .Bc«ll-Bid)th«B 

■1  (D^Ci-Bcell-t-aid) 
•Bd  if 

102  COBtiBB* 

call  v«c_$icopy(Bl,B,BC«ll) 
106  COBtiBB* 

do  500  ir«l,Brif 
do  110  i«l,Bcol 
kl(i)»«(i) 

110  COBtiBBO 

do  120  i«l,BC«2 
k2(i)«B(Be«l-t'i) 

120  COBtiBBO 

do  130  i«BCo2,l,-l 
■(2«i-l)«k2(i) 

130  COBtiBBO 

do  140  isBCol.l.-l 
■(2*i)=kl(i) 

140  COBtiBBO 

600  COBtiBBO 

100  COBtiBBO 
rotBTB 
OBd 


c— —iBitializatioB  of  a  roal  2D  array 

SBbroBtiBo  riBit2(iroB, icol ,a, valBo) 

roal  a(iro«, icol) .tbIbo 

do  5  j=l,icol 
do  10  i=l,iroB 
a(i. j)=valBo 
10  COBtiBBO 

6  COBtiBBO 

rotBTB 

OBd 

c - iaitializatioB  of  a  roal  3D  array 

BBbroBtiao  riBit3 (irov , icol ,lay , a, valBo) 

roal  a( irov, icol, lay), valBo 

do  15  l«l,lay 
do  5  jsl,icol 


do  10  i*l.irow 

10  coat inn* 

S  continno 

15  continno 
rotnrn 
oad 

c — initialization  of  an  intogor  20  array 

snbrout ino  init2 ( iron , icol , ia , valno) 

integer  iaCiron, icol), value 

do  5  jal.lcol 
do  10  isl.iron 
ia(i, j)svalne 
10  continue 

S  continue 
return 
and 

c — -copy  Iron  one  integer  20  array  to  another 

subrout ino  natz.copy ( iron , icol , iu , iv ) 

integer  iu( iron , icol) , iv (iron, icol) 

do  6  jal.icol 
do  10  isl.irov 
iv(i,j)=iu(i,j) 

10  continue 

5  continue 

return 
end 

c — copy  iron  one  real  20  array  to  another 

subroutine  matz.copy(irov,icol,u,v) 

real  u(irov, icol)  ,v(iz'ov, icol) 

do  5  jsl,icol 
do  10  isl.irov 

10  continue 

5  continue 

return 
end 

c - copy  Iron  one  real  30  array  to  another 


■ubront ina  x«atx3_copy ( iro* , icol , lay , v) 
real  a ( iro« . icol , lay) , v( irov . icol , lay) 
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do  IS  iBl.lay 
do  5  jsl.icol 
do  10  isl.irow 

v(i,j,l)=u(i,j,l) 

10  continao 

5  continao 

15  continao 
rotarn 
and 

c - sabroatina  to  choosa  radii  randoaly  lor  poly-diaparaa  syatoa 

sabroatina  r.assign(nl,n2,n3,nl0,n20.n30,rl,r2,r3,pl,p2.iaaad,r) 

implicit  real  (a-b,o-z) 
iaplicit  intogar  (i-n) 

if  (Cm.  It.  nlO)  .and.  (n2.lt  .n20)  .and.  (n3.1t.n30))than 
ra«ran(iaaad) 
if (ra . It . pi ) than 
r=rl 
nl®nl+l 

also  if ((ra.gt.pl). and. (ra.lt. p2))than 
r«r2 
n2»n2't-l 
also 
r»r3 
n3^3-i-l 
end  if 

also  if ( (n2 . It . n20) . and . (n3 . It .n30) ) than 
rasranCiaaad) 

if (ra.lt. (p2-pl)/(1.0-pl))than 
r®r2 
n2=n2+l 
also 
r®r3 
n3=n3+l 
and  if 

alaa  if ((nl.lt. nlO). and. (n2.lt .n20))than 
rasran(iaead) 
if (ra.lt . pl/p2)than 
r*rl 
nl=nl-t’l 
also 
r=r2 
n2^2-*-l 
and  if 

alsa  if ( (nl . It . nlO) . and. (n3 . It .n30))than 
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ra«rui(is««d) 

il  (ra.  It  .pl/(l .  0-p2'«'pl))tb«n 
rarl 
nl*nl+l 

•la* 

r»r3 
n3sn3-t-l 
•nd  H 

•la* 

if (n3 . It . nSOlthen 
r=r3 
ii3=ii3+l 

•la*  il(n2.1t.n20)th*n 
r=r2 
n2=n2+l 

•la* 

r=rl 

•nd  if 
•nd  if 

return 

•nd 


c - calculating  atraaa  tansor 

aubrout in«  atr esa (dm , npt , nonop , epaa , nx , ny , nz , Inl , f t 1 , r , 
Tolun*,a,aigna,iOorl,iaigB) 

inplicit  r«al  (a-h,o-z) 
inplicit  integer  (i-n) 

integer  dn 

real  nx(npt,npt) ,ny(npt,npt),nz{npt,npt) ,lal(npt,npt) 
real  nxx , nyy , nzz , r (npt ) . f t 1 (npt , npt , dm) , a (dn , dm) 

do  877  i=l,npt-l 
do  876  j=i+l,npt 
nx(j,i)=-nx(i,j) 
ny(j,i)=-ny(i,j) 
nz(j,i)=-nz(i,j) 
fnl(j,i)=ynl(i,j) 
do  874  k=:l,dB 

fti( j ,i,k)=-ltl(i, j ,k) 

874  continue 

875  continue 

877  continue 

do  10  isl.dn 
do  15  j=l,dm 
■(jti)=0*0 
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15  continue 

10  continue 

il (nonop. eq. l)then 
do  121  npsl.npt-l 
do  122  kn=np-t-l  ,npt 

il(lnl(np.kn).ge.-ep8a)go  to  122 
nxzsiuc(np,kn) 
nyy®ny(np,kn) 
nzz=nz(np,kn) 
f nlx-f nl (np , kn) *nxz 
Inlyslnl (np , kn) *nyy 
f nlz=f nl (np , kn) *nzz 

8 ( 1 , 1 ) =8 ( 1 , 1 )+r (np) *nxx* (f nlx+f t 1 (np , kn , 1 ) ) 

8 ( 2 , 2 ) =8 (2 , 2) tr (np) onyy * (Inly +f t 1 (np , kn , 2) ) 

8  (  3 , 3) =8  (3 , 3) -t-r  (np)  *nzz*  (f  nlz-f  1 1 1  (np .  kn ,  3  )  ) 

if (iOorl . eq. l)then 

8(1.2)=8(1 , 2)+r (iq>) *nxx* (f nly+ltl (np ,kn . 2) ) 
8 ( 1 , 3) =8 ( 1 , 3)+r (np) *nxx* (Inlz+f t 1 (np ,kn , 3) ) 
8  (2 . 1 )  =8  (2 , 1  )-i-r  (np)  *nyy8  (f  nlx-t-f  tl  (np  ,kn ,  1 )  ) 
8(2,3)=8(2,3)-«'r(iqp)*nyy*(fnlz-t>ftl(np,kn,3)) 
8(3,1)=8(3,  D-i-r  (np)  *nzz*  (f  nlx-*-f  tl  (np  ,kn ,  1 )  ) 
8  (3, 2)=8  (3 .2)-*'r  (iq>)  anzz*  (f  nlyf  f  tl  (np  ,kn ,  2)) 

end  il 

122  continue 

121  continue 

else 

do  721  np=l,npt 
do  722  kn=l,npt 

if (lnl(np,kn) .ge.-ep8a.or.np.eq.kn)go  to  722 
nxx=nx(np,kn) 
nyy=ny(np,kn) 
nzz=nz(np,kn) 
f nlx=f nl (np , kn) *nxx 
fnly=fnl(np,kn)*nyy 
f nlz=f nl (np , kn) enzz 

8(l,l)=8(l,l)+r(np)*nxx*(lnlx+ftl(np,kn,l)) 

8(2,2)=8(2,2)+r(np)*nyy*(fnly+ltl(np,kn,2)) 

8  (3 , 3)=8  (3 , 3)'t’r  (np)*nzz*  (f  nlz-ff  1 1  (np  ,kn ,  3)  ) 
if (iOorl . eq. l)then 

8 ( 1 ,2)=8(1 ,2)+r(np)*nxx* (fnly+f tl (np ,kn , 2)) 
8(1,3)=8(1, 3)-i-r (np)enzze (f nlz'ff tl  (np ,kn , 3) ) 
8(2,l)=8(2,l)+r(np)*nyy*(fnlx+ftl(np,kn,l)) 
8(2,3)=8(2,3)-*-r  (np)  *nyy*  (f  nlz+f  tl  (np  ,kn ,  3)  ) 
s(3,  1)=8(3,  l)4r(np)enzze(fnlx-t'ftl(np,kn,  1)) 
8(3, 2) =8 (3 , 2)+r (np) *nzz* (f nly+f tl (np ,kn , 2) ) 


end  if 


continns 

continu* 

•nd  if 

do  20  isl.dM 
do  25  jsl.da 

il (monop. oq. l)s ( j , 1)32. 0*t(j , i)/voluao 
il(moBop.oq.2)s(j,i)38(j.i)/voliiBO 
continue 
continue 

if (ie ign . eq . 2) then 

•ign&3-(>(l,l)4a(2.2))/2.0 

else 

sigmn3-(s(l,l)-«-s(2.2)-t-s(3,3))*0.333333 
end  if 

return 

end 


subroutine  conduct(npt,Bznc,nnb,dx,dy,dz,sBtm,x,y,z,r, 
nx , ny , nz , cond , brx , bry , brz , cur , voltl , c_e« , 
cm ,  eres ,  cdm  .nap  .nappairt  ,adjc  ,kxx  ,kyy ,  kzz ,  grad ,  berk) 

real  sstm(3,3) ,x(npt) ,y(npt),z(iq>t) ,r(npt) 
real  nx(npt,npt) ,By(npt,npt),nz(npt,npt) 
real  cond(npt ,npt) ,brx(npt) ,bry(npt) ,brz(npt) 
real  cur (npt , npt ) , volt 1 (npt ) , c_e« (3,3) 
real  cm(npt)  .cresCnpt)  ,cdm(npt) 

integer  Bap(npt) ,kxx(mmc,nnb) ,kyy(mxnc,nnb) ,kzz(mxnc,nnb) 
integer  adjc(nxnc,nnb) ,nappart(Bmc) 

call  rinit2(3,3,c_ell,0.0) 

call  relax_v(npt , cond,brx, cm, cres , cdm) 
call  vec_$copy(cm,Toltf  ,npt) 

Brite(60,*)ir,ist, ’  x’ 

lx=0.0 

fy=0.0 

fz^O.O 

do  41  np=l,npt 
snm=0.0 

do  42  j_kn=l,nnb 

kn3Bappart(adjc(Bap(np) , j_kn)) 
if (kn.eq.O)go  to  42 


cur(np,kn)30.0 
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il(abt(cond(ap,ka)).lt.l.0«-16.or.&p.«q.kB)go  to  42 

xl«)uz(Mip(Bp) .  j_]m)*dZ't'k7j(Bap(&p) .  j_kn)*tan(MtrB(2,  l))*dy 
■•■kzz  (b^  (np) ,  j  _kn)  *tui(satza(3 . 1 )  )  *dz 

zkn«z(kii}-zl 

cur  (up ,  kn)  Bcoad  (np ,  kn)  *  (  (x(Bp)-zka)  *barR*grad 
♦toltf (np) -Toltl (kn) ) 
nrito (58 , *) cor (np,kn) 
sna«siia+cnr(np.kn) 

tx*ix*x (np) *bnrR*nx (np , kn) *cnr (np , kn) 
ijsi  y-tr  (np)  *barR*ny  (np ,  kn)  *cnr  (np .  kn) 
f  zef  x-fr  (np)  *bnrR*nz  (np .  kn)  *ciir  (np .  kn) 
continno 

nrit«(60 ,  *)iq> ,  aim 
continue 

ix^x/  (dz*dy*dz*burR*barR*barR) 

tj^j/  (dz*dy*dz*barR*barR*bnrR) 

lz^z/(dx*dy*dz*barR*barR*barR) 

c_a«(l,l)»-fx 

c_e«f(1.2)*-fy 

C_6fl(l,3)a-f2 

call  relax. V (npt , cond , bry , cxn, cr es , cdxn) 
call  vec_$copy(cxn.voltl,npt) 

«rite(60,*)ir,ist, ’  y‘ 

ix*0.0 

ly=0.0 

Xz=0.0 

do  51  np-l,npt 
auBsO . 0 

do  52  j_knsl,nnb 

knsmappart(adjc(Bap(np) , j_kn)) 
iX(kn.eq.O)go  to  52 

cur(np,kn}=0.0 

iX(abs(cond(np,kn)) .lt.l.0e-15.or.np.eq.kn)go  to  52 

ylskyy(nap(np) ,  j  Jcn)*dy*'kxx(i>ap(np) ,  j_kn)*tan(sstm(l  ,2))*dx 
■^kzz  (nap(np) ,  j_kn)  *tan(88tm(3 , 2)  )  *dz 

ykn=y(kn)-yl 

cur (np , kn) «cond (np , kn) * ( (y (np) -ykn) ^barRegrad 
+voltX (np) -voltt (kn) ) 
suBsaun+cur (np , kn) 

Ixsfx+r (np) *barR*nx (np , kn) *cur (np , kn) 
f y=f y+r (np) *barR*ny (np, kn) ♦cur (np , kn) 

IzsXz-t'r  (np)  ♦barRenz(np,kn)  ecur  (np ,  kn) 
continue 


*r it • ( 60 , * } up , ana 
continna 

ix^x/  (dz*dy*dz*barR«barR*burR) 

(dx*dy*dzabarR*barR*barR) 
f  z^z/  (dz*dy*(tz*barR*baxR*barR) 
c_a«(2,l)*-lx 
c.aff(2,2)s-fy 
c_e«(2,3)=-fz 

call  r claz.v (apt , cond , brz , czn , cr as , cdzn) 
call  Tac_$copy(cxn.voltl,npt) 

«rita(60,a)ir,iat. *  z* 

lx»0.0 

fy=0.0 

lz=0.0 

c_top=0.0 

c_bot=0.0 

do  61  npsl.npt 
suasO . 0 

do  62  j_kn=l,tmb 

kn=aappart(adjc(aap(iip) ,  j_kn)) 
il(kn.aq.0)go  to  62 

cur(xip,kn)s0.0 

il(abs(cond(np,kii)).lt.l.Oa-lS.or.np.eq.kn)go  to  62 

zlskzzCaapCnp) , j_kn}adz-^kxx(Bap(np) , j_lm)atan(88trn(l ,3))adx 
+kyy  (map  (np) ,  j  _kii)*tan  (s8tzii(2 , 3)  )  *dy 

zk&=z(kii}-zl 

cur (np , kn) =cond (np , kn)* ( (z(np}-zkn) *barR*grad 
+voltl (np) -voltl (kn) ) 


it(kzz(nap(np) ,nap(kn)) .aq.l)c_top=c_top+cur(np,kn) 
il(kzz(nap(np) ,Bap(kn)) .aq.-l)c_bot=c_bot+cur(np,kn) 


il(np.la.l0)n:ita(60,*)np,kn,cur(np,kn) 
suassua-fcur  (np ,  kn) 

f x=lx+r (np) abarR*nx (np ,kn) *cur (np , kn) 

1 y =ly +r (np) abarRany (np, kn) acur (np , kn) 

1 z=f zar (np) abarRanz (np, kn) acnr (np , kn) 
continue 

nrita(60 , a)np , 8ua 

«rite(69 , a)z(np) ,z(np)agradpavoltt (np) 


continue 
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ix*tx/  (dz*d]r*dz*barR*barR«biLrll) 
(dz*dy*dz*bu:R*barR*barR) 
ix*tz/ (dz*dj*4z*barR*barR*barR) 
c_«ll(3,l)»-lx 
c.««(3.2)*-fy 
c_««(3,3)— 1* 

c - 

c  «rit«(96,*)c_top,c_bot 

c  «r it • ( 95 , * ) dz*barR*  e_t op/ (dx*bar R»dy *barR*grad*dz*barR) 

c  «rit«(9S,*)c_«ll(3,3) 

return 

end 


c - linear  equation  aolver 

subroutine  relax_v(npt ,cond,br,cxn,cres,cdxn) 
inplicit  real  (a-h.o-z) 
iaplicit  integer  (i-n) 

real  cond(npt,npt) .br(npt) ,cxn(npt) ,cros(npt) ,cdxn(npt) 
xerosQ . 0 

call  vec.$init(czn,npt,zero) 
cadi  vec.$init(cdxn,npt,zero) 
call  vec_$init(cre8,npt,zaro) 

inuB=0 

100  inun=inun-fl 

raiax=0.0 
inax=0 
cresun^O.O 
if (inun . gt . 1 }then 
do  22  i=l,npt 

il(ab8(cond(i, jnax)).lt.l.0e-15)go  to  21 
cr  as  (  i  )  seres  (  i  ) -fcondC  i ,  jnax)  *cdxn  (jnaz  ) 

21  cresuBscresun-fabs(cres(i)) 
if (abs(cres(i)) .gt.raaxlthen 

naxsabs  (  cr  es  (  i  )  } 
inaxsi 
end  if 

22  continue 

else 

do  20  isl.npt 

if (br(i).eq.0.0)go  to  20 
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crM(i)«-br(i) 
cTMiui*crMiia-»^ab«  (cr«s(i)  ) 
if (abt (  arts ( i) ). gt . xmax) than 
nukzsaba  (  cr  •>  (  i  )  ) 
imaxsi 
•nd  if 

20  continue 

crasual^crasna 
and  if 

jMXsinax 

cdzn  (  j  BU  )  »-cr  as  (  iaax) /cond(  inax ,  jaax) 
cxn(jBax)scxn(jBax)+cdxn(jBax) 

if (crasna.gt. crasuBlaO. 00001. and. inuB.lt. 10000)go  to  100 
c  vrita(60,*) ’crasuBl.crasuB.inuB* ,crosuBl.crasuB,inuB 

return 

and 


c - copy  one  half  of  a  20  array(aboTa  diagnol)  to  another  half 

subroutine  copyhalf (irou.u) 

real  u(irov,irou) 

do  6  i=l,irow-l 
do  10  j=i+l,irow 
u(j,i)=u(i,j) 

10  continue 
S  continue 
return 
end 


c — 


subroutine  degree_phi(nz,kphi) 
inplicit  real  (a-b,o-z) 
inplicit  integer  (i-n) 

real  nz 

integer  kphi(12) 

dnz=0. 1666667 
if (nz .gt . 0 . 0)then 
k=6-int(nz/dnz) 

else 

k=-int  (nz/dnz)  7 
end  if 

kphi(k)  skphi  (k) 
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return 

•nd 


isbrovtiiia  d«{r««_tli«ta(iuiz,iui]r.kth«ta) 
iaplicit  r«ml  (a-h.o-x) 
iaplicit  intagar  (i-a) 

raal  u.&j.niuc.imy 
intagar  ktliata(24) 

pis3. 141693 

nx*nnz/aqrt  (muaanz-t^imjaiuij) 
aysnny/aqrt  (anxaaax't'aByaBBy) 

if  (ay .  ga .  0 . 0)tliaa 
if (ax . ga . 0 . 0)thaa 

thataBaaia(aba (ay) ) 

alaa 

thataspi‘-asia(abs(ay)  ) 
aad  if 

alaa 

if (ax . ga . 0 . 0)thaa 

tliatas2 .  Oapl-aala(abs  (ay)  ) 

alaa 

thataspi'*'aaia(aba  (ay)  ) 
aad  if 
aad  if 

thataa(tbata/pi)*180.0 
kxiat  (tbata/15 . 0)'f  1 
kthata(k)-kthata(k)+l 

ratara 

aad 
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